{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multi-point stress approximation (MPSA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Porepy supports mpsa discretization for linear elasticity problem:\n",
    "\\begin{equation}\n",
    "\\nabla\\cdot \\sigma = -\\vec f,\\quad \\vec x \\in \\Omega\n",
    "\\end{equation}\n",
    "where $\\vec f$ is a body force, and the stress $\\sigma$ is given as a linear function of the displacement\n",
    "\\begin{equation}\n",
    "\\sigma = C:\\vec u.\n",
    "\\end{equation}\n",
    "\n",
    "The convention in porepy is that tension is positive. This means that the Cartesian component of the traction $\\vec T = \\sigma \\cdot \\vec n$, for a direction $\\vec r$ is positive number if the inner product $\\vec T\\cdot \\vec r$ is positive. The displacements will give the difference between the initial state of the rock and the deformed state. If we consider a point in its initial state $\\vec x \\in \\Omega$ and let $\\vec x^* \\in \\Omega$ be the same point in the deformed state, to be consistent with the convention we used for traction, the displacements are given by $\\vec u = \\vec x^* - \\vec x$, that is, $u$ points from the initial state to the finial state.\n",
    "\n",
    "To close the system we also need to define a set of boundary conditions. Here we have three posibilities, Neumman conditions, Dirichlet conditions or Robin conditions, and we divide the boundary into three disjont sets $\\Gamma_N$, $\\Gamma_D$ and $\\Gamma_R$ for the three different types of boundary conditions\n",
    "\\begin{equation}\n",
    "\\vec u = g_D, \\quad \\vec x \\in \\Gamma_D\\\\\n",
    "\\sigma\\cdot n = g_N, \\quad \\vec x \\in \\Gamma_N\\\\\n",
    "\\sigma\\cdot n + W \\vec u= g_R,\\quad \\vec x\\in \\Gamma_R\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "To solve this system we first have to create the grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import porepy as pp\n",
    "\n",
    "# Create grid\n",
    "n = 5\n",
    "g = pp.CartGrid([n,n])\n",
    "g.compute_geometry()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to define the stress tensor $C$. In porepy the constutitive law,\n",
    "\\begin{equation}\n",
    "\\sigma = C:u = 2  \\mu  \\epsilon +\\lambda  \\text{trace}(\\epsilon) I, \\quad \\epsilon = \\frac{1}{2}(\\nabla u + (\\nabla u)^\\top)\n",
    "\\end{equation}\n",
    "is implemented, and to get the tensor for this law we call:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create stiffness matrix\n",
    "lam = np.ones(g.num_cells)\n",
    "mu = np.ones(g.num_cells)\n",
    "C = pp.FourthOrderTensor(mu, lam)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we need to define boundary conditions. We set the bottom boundary as a Dirichlet boundary, and the other boundaries are set to Neuman."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define boundary type\n",
    "dirich = np.ravel(np.argwhere(g.face_centers[1] < 1e-10))\n",
    "bound = pp.BoundaryConditionVectorial(g, dirich, ['dir']*dirich.size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We discretize the stresses by using the multi-point stress approximation (for details, please see: E. Keilegavlen and J. M. Nordbotten. “Finite volume methods for elasticity with weak symmetry”. In: International Journal for Numerical Methods in Engineering (2017))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now define the values we put on the boundaries. We clamp the bottom boundary, and push down by a constant force on the top boundary. Note that the value of the Neumann condition given on a face $\\pi$ is the integrated traction $\\int_\\pi g_N d\\vec x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_faces = np.ravel(np.argwhere(g.face_centers[1] > n - 1e-10))\n",
    "bot_faces = np.ravel(np.argwhere(g.face_centers[1] < 1e-10))\n",
    "\n",
    "u_b = np.zeros((g.dim, g.num_faces))\n",
    "u_b[1, top_faces] = -1 * g.face_areas[top_faces]\n",
    "u_b[:, bot_faces] = 0\n",
    "\n",
    "u_b = u_b.ravel('F')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We discretize this system using the Mpsa class. We assume zero body forces $f=0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "parameter_keyword = \"mechanics\"\n",
    "\n",
    "mpsa_class = pp.Mpsa(parameter_keyword)\n",
    "f = np.zeros(g.dim * g.num_cells)\n",
    "\n",
    "specified_parameters = {\"fourth_order_tensor\": C, \"source\": f, \"bc\": bound, \"bc_values\": u_b}\n",
    "data = pp.initialize_default_data(g, {}, parameter_keyword, specified_parameters)\n",
    "mpsa_class.discretize(g, data)\n",
    "A, b = mpsa_class.assemble_matrix_rhs(g, data)\n",
    "\n",
    "u = np.linalg.solve(A.A, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can plott the y_displacement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/eke001/python_env/eris/lib/python3.6/site-packages/mpl_toolkits/mplot3d/axes3d.py:753: UserWarning: Attempting to set identical bottom==top results\n",
      "in singular transformations; automatically expanding.\n",
      "bottom=0.0, top=0.0\n",
      "  'bottom=%s, top=%s') % (bottom, top))\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxUAAAKhCAYAAAAmMYPdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xl4VeW99vH7yUAgA2MAQxIZCsikIAkoaq0iDrUWi1LFSoWDWHtKr6NH3x4HlOJw1NpatQfb89rqcaBO7VvRgxhlsINADEGIQJRJA4Qwy5R5D8/7R5JdkuwkOyx21gr5fq5rX+6119pr//YyxvXL/TxrGWutAAAAAOBkxbhdAAAAAID2jaYCAAAAgCM0FQAAAAAcoakAAAAA4AhNBQAAAABHaCoAAAAAOEJTAQAAAMARmgoAAAAAjtBUAAAAAHAkroX13G4bAAAAbcG4XUAkBhtjy138/D3SB9baq1wsIayWmgoAAAAAtcol3e7i58+XUl38+CYx/AkAAACAIyQVAAAAQISMOIEOh6QCAAAAgCM0FQAAAAAcIb0BAAAAImQkxbtdhAeRVAAAAABwhKQCAAAAiBATtcMjqQAAAADgCE0FAAAAAEdIbwAAAIAIMVE7PJIKAAAAoAMwxvQ0xiw1xmyt/WePMNuMMcasNsZsMsZ8Zoy5MZJ901QAAAAAEaqbqO3Ww6F7JS231g6RtLx2uaFySbdYa0dKukrSM8aY7i3tmKYCAAAA6BiulfRy7fOXJX2v4QbW2i3W2q21z0sk7ZfUu6UdM6cCAAAAaD9SjTH5Jyw/b619PsL39rXW7ql9vldS3+Y2NsaMl9RJ0vaWdkxTAQAAAETIAxO1D1prs5taaYxZJumMMKvmnrhgrbXGGNvMftIkvSpphrU22FJRNBUAAADAacJaO6mpdcaYfcaYNGvtntqmYX8T23WV9J6kudba3Eg+lzkVAAAAQITa+UTtdyXNqH0+Q9I7jb6fMZ0kvS3pFWvtnyPdMU0FAAAA0DE8IelyY8xWSZNql2WMyTbG/KF2mxskXSxppjFmfe1jTEs7NtY2OZRKkppdCQAAAJwixu0CIjHQGPtzFz//X6S1zc2pcAtzKgAAAIAIeWCiticx/AkAAACAIyQVAAAAQITqJmqjPpIKAAAAAI7QVAAAAABwhPQGAAAAiBATtcMjqQAAAADgCE0FAAAAAEcY/gQAAABEiOFP4ZFUAAAAAHCEpAIAAABoBU6gGyOpAAAAAOAITQUAAAAAR0hvAAAAgAgxUTs8kgoAAAAAjpBUAAAAABEy4gQ6HJIKAJ6WnJysL7/8Muy6l156SRdddFEbVwQAABqiqQCg1157TdnZ2UpOTlZaWpq+/e1v6+OPPz7p/RljtG3btnqvHT9+XHfddZcGDBigpKQknXnmmZo6dao++eSTZvdVWlqqQYMGnVQd1dXVmj9/voYMGaKkpCQNGDBAs2bNUlFR0UntL1rmz5+v6dOnu11GVBUVFckYI7/f73YpAIAooKkAOrhf//rXuvPOO3X//fdr37592rlzp37yk5/onXfeafW+mjphrKqq0sSJE7VhwwYtXrxYx44d0+eff65p06bp/fffb9W+WmPq1Kl699139dprr+no0aMqKChQVlaWli9f7njfOPVoOAC0B3UTtd16eBVNBdCBHT16VPPmzdNzzz2n6667TklJSYqPj9d3v/td/fKXv5Qk5eXlacKECerevbvS0tL005/+VNXV1aF9GGP03HPPaciQIRoyZIguvvhiSdLo0aOVnJysN998U6+++qqKi4u1aNEijRo1SrGxsUpKStLUqVM1f/78JvdV91pd6nHo0CFNnjxZXbt21fjx47V9+/Ymv9uyZcu0dOlSvfPOOxo3bpzi4uLUrVs3zZkzR7feeqskqaSkRJMnT1bPnj01ePBg/f73vw+9f/78+fr+97+v6dOnKyUlRWeffba2bNmixx9/XH369FFmZqY+/PDD0PaXXHKJ7rvvPo0fP15du3bVtddeq6+//lqS9Ne//lUZGRn16hswYICWLVumnJwcPfbYY3rzzTeVnJys0aNHh/7d3HrrrUpLS1N6eroeeOABBQKBsN+1oqJCM2bMUI8ePTR8+HA9+eST9T6vpKRE119/vXr37q2BAwfqN7/5TWhdVVWV7rzzTvXr10/9+vXTnXfeqaqqqnp1P/nkk+rTp4/S0tK0aNEiLVmyREOHDlXPnj312GOPhfYVDAb1xBNP6Bvf+IZ69eqlG264IXQM6n4uunfvruTkZK1evVovvfSSLrzwQv37v/+7evXqpXnz5qlnz57asGFDaJ/79+9XYmKiDhw40OS/awCA+2gqgA5s9erVqqys1JQpU5rcJjY2Vk8//bQOHjyo1atXa/ny5frtb39bb5tFixbpk08+UWFhof7+979LkgoKClRaWqobb7xRy5Yt05VXXqmkpKQWazpxXw3NmTNHnTt31p49e/Tiiy/qxRdfbHI/y5Yt0/jx45WZmdnkNtOmTVNGRoZKSkr05z//Wffff79WrFgRWv+///u/+uEPf6jDhw/r3HPP1ZVXXqlgMKjdu3dr3rx5uv322+vt75VXXtGLL76oPXv2KC4uTv/2b//W4ve96qqrdP/99+vGG29UaWmpCgoKJEkzZ85UXFyctm3bpnXr1unDDz/UH/7wh7D7eOihh1RUVKQvv/xSS5cu1cKFC0PrgsGgvvvd72r06NHavXu3li9frmeeeUYffPCBJOk///M/lZubq/Xr16ugoEB5eXl69NFHQ+/fu3evKisrtXv3bj388MO67bbbtHDhQq1du1b/+Mc/9Mgjj+irr76SJP3Xf/2XFi1apL/97W8qKSlRjx49NGfOHEkK/VwcOXJEpaWlmjBhgiTpk08+0aBBg7Rv3z49+OCDmjZtWr36X3/9dV122WXq3bt3i8cSANpC3URttx6eZa1t7gHgNLZw4ULbt2/fVr3n6aeftt/73vdCy5Ls8uXL620jyW7dujW0fNlll9l77rkntLxu3TrbrVs3m5KSYocOHRrRvvx+v42Li7Off/55aN19991nL7zwwrB1zp492954441Nfo+dO3famJgYe+zYsdBr9957r50xY4a11tqf//zndtKkSaF17777rk1KSrJ+v99aa+2xY8esJHv48GFrrbXf+ta36n3HTZs22fj4eOv3++1HH31k09PT631+//797dKlS0OfdfPNN4fW7d2713bq1MmWl5eHXnvttdfsJZdcEva7DBw40Obk5ISWf//734c+Lzc312ZmZtbb/rHHHrMzZ8601lo7aNAg+95774XW5eTk2P79+1trrf3oo49s586dG33n3Nzc0PZjx461b7/9trXW2mHDhtlly5aF1pWUlNi4uDjr8/nsV199ZSVZn88XWv8///M/jWqrqzcYDFprrc3KyrJvvvlm2O8N4LTT0nmpJx7DJJvr4kNSvtvHINzD0w0PgOjq1auXDh48KL/fr7i48L8OtmzZorvuukv5+fkqLy+X3+9XVlZWvW2aSwPqPmfPnj2h5TFjxujIkSNatmyZZs+eHdG+Dhw4IL/fX299//79m/3MLVu2NLm+pKREPXv2VEpKSr395efnh5b79u0bet6lSxelpqYqNjY2tCzVTCTv3r17o9r79+8vn8+ngwcPNllDU3bs2CGfz6e0tLTQa8FgsMljU1JSUm/dic937NihkpKSUI2SFAgE9M1vfjP03hOPY//+/VVSUhJa7tWrV6Pv3PC4lJaWhj5rypQpion5ZwgeGxurffv2NfldG36n8847T4mJifrrX/+qtLQ0bdu2TZMnT27y/QAAb2D4E9CBTZgwQQkJCVq0aFGT2/zrv/6rhg0bpq1bt+rYsWN67LHHZK2tt40xptnPueyyy/Thhx+qrKysxZqa2lfv3r0VFxenXbt2hV7buXNnk/uZNGmS8vLyVFxcHHZ9v3799PXXX+v48eP19peent5ijU1pWFt8fLxSU1OVlJSk8vLy0LpAIFBvjkDD75yZmamEhAQdPHhQR44c0ZEjR3Ts2DFt2rQp7OempaXV+54n1pGZmamBAweG9nPkyBEdP35cS5YskVRzHHbs2FGv7n79+p3U98/MzNT7779f77MqKyuVnp7e5L/XcK/PmDFDCxcu1KuvvqqpU6eqc+fOJ1UPAEQDE7XDo6kAOrBu3brp4Ycf1pw5c7Ro0SKVl5fL5/Pp/fff13/8x39IqrkUbNeuXZWcnKwvvvhCv/vd71rcb9++fevdW+KWW25RWlqapkyZoo0bNyoQCKiysrJeKtCS2NhYXXfddZo/f77Ky8tVWFiol19+ucntJ02apMsvv1xTpkzR2rVr5ff7dfz4cf33f/+3XnzxRWVmZuqCCy7Qfffdp8rKSn322Wd64YUXHF3adeHChSosLFR5ebnmzZunqVOnKjY2VkOHDlVlZaXee+89+Xw+Pfroo6HJ0FLN8SoqKlIwGJRU0yRcccUVuvvuu3Xs2DEFg0Ft375df/vb38J+7g033KDHH39chw8f1u7du7VgwYLQuvHjxyslJUW/+MUvVFFRoUAgoI0bN2rNmjWSpJtuukmPPvqoDhw4oIMHD+rhhx8+6WPw4x//WHPnzg01KQcOHAhdRax3796KiYlp8p4jJ5o+fbrefvttLVy4ULfccstJ1QIAaFs0FUAHd/fdd+vXv/61Hn30UfXu3VuZmZlasGCBvve970mSfvWrX+m1115TSkqKbrvtNt14440t7nP+/PmaMWOGunfvrrfeekudO3fWRx99pBEjRug73/mOunbtqrPOOktr1qzRW2+9FXGtCxYsUGlpqc444wzNnDlT//Iv/9Ls9n/+85919dVX68Ybb1S3bt00atQo5efna9KkSZJqJgEXFRWpX79+mjJlih566KHQupPxwx/+UDNnztQZZ5yhysrK0FWWunXrpt/+9reaPXu20tPTlZSUVO/qTN///vcl1Qw1Gjt2rKSaSd/V1dUaMWKEevTooalTp9YbQnaiefPmKSMjQwMHDtSkSZM0depUJSQkSKppxhYvXqz169dr4MCBSk1N1ezZs3X06FFJ0gMPPKDs7Gydc845OvvsszV27Fg98MADJ/X977jjDk2ePFlXXHGFUlJSdP7554fuQ5KYmKi5c+fqwgsvVPfu3ZWbm9vkfjIzMzV27FgZY0LDtADAK5ioHZ5pOIyhgWZXAgBqXHLJJZo+fXqjOSJu+N3vfqc33nijyWSjPZg1a5b69etX70pUAE57zY+l9YiRxtjXXfz80dJaa222iyWERVIBAO3cnj17tHLlSgWDQW3evFlPPfVUs5cJ9rqioiL95S9/Cd1PBADgfTQVANDOVVdX6/bbb1dKSoomTpyoa6+9Vj/5yU/cLuukPPjggxo1apR+9rOfaeDAgW6XAwCNMFE7PIY/AQAAwAvaxfCnUcbYP7n4+SM8OvzJy/M9AAAAAE+pm6iN+hj+BAAAAMARmgoAAAAAjpDeAAAAABGqm6iN+kgqAAAAADhCUwEAAADAEYY/AQAAABFi+FN4JBUAAAAAHCGpAAAAAFqBE+jGSCoAAAAAOEJTAQAAAMAR0hsAAAAgQkZSvJtn0H4XP7sZJBUAAAAAHCGpAAAAACJkjBRHUtEISQUAAAAAR2gqAAAAADjC8CcAAAAgQsZI8bFuV+E9JBUAAAAAHCGpAAAAACLk+kRtjyKpAAAAAOAITQUAAAAARwhvAAAAgAi5fkdtjyKpAAAAAOAIfRYAAAAQKSOJS8o2QlIBAAAAwBGaCgAAAACOMPwJAAAAiJQRZ9BhkFQAAAAAcIQ+CwAAAIgUSUVYJBUAAAAAHKGpAAAAAOAI4Q0AAADQGpxBN0JSAQAAAMARmgoAAAAAjhDeAAAAAJEykmLdLsJ7SCoAAAAAOEJSAQAAAESK+1SERVIBAAAAwBGaCgAAAACOEN4AAAAAkWL4U1gkFQAAAAAcoc8CAAAAWoNLyjZCUgEAAADAEZoKAAAAAI4w/AkAAACIFBO1wyKpAAAAAOAIfRYAAAAQKZKKsEgqAAAAgA7AGNPTGLPUGLO19p89mtm2qzGm2BizIJJ901QAAAAAHcO9kpZba4dIWl673JRHJP090h3TVAAAAACtEeviw5lrJb1c+/xlSd8Lt5ExJktSX0kfRrpjmgoAAACgY+hrrd1T+3yvahqHeowxMZKekvR/WrNjppkAAAAAkXJ/onaqMSb/hOXnrbXP1y0YY5ZJOiPM++aeuGCttcYYG2a7n0haYq0tNsZEXBRNBQAAANB+HLTWZje10lo7qal1xph9xpg0a+0eY0yapP1hNpsg6ZvGmJ9ISpbUyRhTaq1tbv4FTQUAAADQQbwraYakJ2r/+U7DDay1N9c9N8bMlJTdUkMh0VQAAAAAkXN/+JMTT0h6yxhzq6Qdkm6QJGNMtqQfW2tnn+yOjbXhhlKFNLsSAAAAOEUiH8DvouyuxuaPd+/zzXKtbW74k1vab58FAAAAtLX2nVREDZeUBQAAAOAITQUAAAAARwhvAAAAgNZwfmfr0w5JBQAAAABHaCoAAAAAOMLwJwAAACBSXP0pLJIKAAAAAI7QZwEAAACRIqkIi6QCAAAAgCM0FQDQRtasWaNzzjlHlZWVKisr08iRI7Vx40a3ywIAwDHCGwBoI+PGjdPkyZP1wAMPqKKiQtOnT9eoUaPcLgsA0BpG3KciDGOtbW59sysBAK1TXV2tcePGqXPnzlq1apViY/k/EwDUMm4XEInsnsbmT3Lv882ftNZam+1eBeGRVABAGzp06JBKS0vl8/lUWVmppKQkt0sCALQGE7XDYk4FALSh22+/XY888ohuvvlm3XPPPW6XAwDAKUGfBQBt5JVXXlF8fLx+8IMfKBAI6IILLtCKFSs0ceJEt0sDAMAR5lQgaoLBoAKBgOLi4mRMuxgmCQAA3NMuThayexmb/x33Pt+86s05FQx/QtRYa+X3+7V9+3a10LwCAACgHWP4E6Ju165dyszMVKdOnUgsAABA+8YlZcMiqUCbCAaDqq6uJrEAAAA4DdFUoE0YY7Rz504aCwAAgNMQw5/QJowxKi4uVkZGhqqrqxkKBQAA2ifuUxEWSQXaVExMjILBIJO3AQAATiP0WWhzMTExTN4GAADtE0lFWCQViBprrQoKChQMBsOuDwaD2rZtG4kFAABAO0dTgagxxujMM89UeXm5ysvLG62PiYlRcXGxqqqqmmw8IpWTk6OzzjpLgwcP1hNPPOFoX6eDWbNmqU+fPho1apTbpXjCrl27dOmll2rEiBEaOXKknn32WbdLcl1lZaXGjx+v0aNHa+TIkfr5z3/udkmeEAgEdO655+qaa65xuxRPGDBggM4++2yNGTNG2dmeu9cWAA+hqUBU9ejRQ126dNG6desUCATCbmOtVXV19Uk3FoFAQHPmzNH777+vwsJCvf766yosLHRSdrs3c+ZM5eTkuF2GZ8TFxempp55SYWGhcnNz9dxzz3X4n5GEhAStWLFCBQUFWr9+vXJycpSbm+t2Wa579tlnNXz4cLfL8JSPPvpI69evV35+vtulAN4R5+LDo2gqEHWxsbEaM2aMKioqdOTIkUbrY2JiZK3V9u3bT6qxyMvL0+DBgzVo0CB16tRJ06ZN0zvvvHMqSm+3Lr74YvXs2dPtMjwjLS1NY8eOlSSlpKRo+PDh2r17t8tVucsYo+TkZEmSz+eTz+fr8PObiouL9d5772n27NlulwIA7Q5NBdpEUlKSEhMTVVhYqIMHDzZaXzcU6mQSi927dyszMzO0nJGR0eFPGNG0oqIirVu3Tuedd57bpbguEAhozJgx6tOnjy6//PIOf0zuvPNOPfnkk4qJ4X+NdYwxuuKKK5SVlaXnn3/e7XIAeBi/OdFmYmJilJWVpa1bt8rn8zW53ckmFkBLSktLdf311+uZZ55R165d3S7HdbGxsVq/fr2Ki4uVl5enjRs3ul2SaxYvXqw+ffooKyvL7VI85eOPP9ann36q999/X88995z+/ve/u10S4D4jKdbFh0fRVKBNJSQkKDs7W9XV1WHThLqb5LVm8nZ6erp27dqlq666SlLNEIb09PRTWnd7tWPHDrdL8Ayfz6frr79eZWVluu6669wuxzOuuuoqde/eXZdeemmHnoezcuVKvfvuuxowYICuvPJKrVixQtOnT3e7LNfV/S695ZZbNGXKFOXl5blcEQCvoqlAm4uPj1diYqL27Nmj6urqsNsYYyJOLMaNG6etW7dq9+7dqq6u1htvvKHJkyef6rLbJb/f73YJnmCt1a233qrhw4eH5hF0dAcOHNCRI0d08OBBVVRUaOnSpRo2bJjbZbnm8ccfV3FxsYqKijRgwABNnDhRCxcudLssV5WVlen48eOSpH379unDDz/kinKA9M/7VDBRux6aCrjCGKOxY8fK7/eHvVdFaxKLuLg4LViwQFu3btXw4cN1ww03aOTIkdEs3/NuuukmTZgwQZWVlcrIyNALL7zgdkmuWrlypV599VWtWLFChYWFGjNmjJYsWeJ2Wa7as2ePLr30UhUWFmrcuHG6/PLLuYwq6tm3b58uuugijR49Wl988YW+853vhBJhAGjIw/0OTncxMTFKTExUZWWlvvjii7Db1CUW3/jGN5qdPHn11Vdr1KhRXPKw1uuvvy5Jys7O5phIuuiii0KNK8ekxjnnnKN169ZxPMJISUnR4sWL3S7DdYMGDVJBQYGkmv9u5s6d63JFALyMpgKuGzlypLZs2aKKiopGqURdYjF27DiVlh5tcV8d95KYsZLC3wek4x6Thv55jDgm9dUcj6Z/hjqi8D8jHfsYtfTfTUpKNx071viy4cBpp274E+rhkMB1xhgNHTpUJSUlob+KNVTTUMxv07ral/mSHnW7CI97QNIv3C7Cw+4RP0MteUD8Hmra8ePz3S4BgItoKuAJxhglJCSoV69e2rJlCxOMAQCAd3n40q5uYaI2POXMM89UfHy81q5d22jyNgAAALyJpgKeEx8fr4EDB6q8vFxVVVVulwMAAIAWMPwJntSnTx8lJCRo7dq13F0bAAB4BxO1wyKpgGfFxcVp5MiRKi8vd7sUAAAANIM+C57WrVs3denSxe0yAAAAapBUhEVSAc+LjeUSCwAAAF5GnwUAAE6J7Ozs0PPU1FTl5OS4WA2AtkRTAQAATon8/Pwm1w0YMEApKSmKjY1VXFxcs9sCnsbwp7A4JAAAoE189NFHSk1NdbsMAFFAUwEAAAC0BtM9G2GiNgAAiDpjjK644gplZWXp+eefd7scAKcYSQUAAIi6jz/+WOnp6dq/f78uv/xyDRs2TBdffLHbZQE4RUgqAABA1KWnp0uS+vTpoylTpigvL8/lioCTVDdR262HR9FUAACAqCorK9Px48dDzz/88EONGjXK5aoAnEoe7ncAAMDpYN++fZoyZYokye/36wc/+IGuuuoql6sCThKXlA2LQwIAAKJq0KBBKigocLsMAFHE8CcAAAAAjpBUAAAAAK3BfSoaIakAAAAA4AhNBQAAAABHGP4EAAAARIqrP4VFUgEAAADAEfosAAAAIFIkFWGRVAAAAABwhKYCAAAAgCOENwAAAECkjLhPRRgkFQAAAAAcIakAAAAAIsVE7bBIKgAAAAA4QlMBAAAAwBHCGwAAAKA1OINuhKQCAAAAgCP0WQAAAECkmKgdFkkFAAAAAEdoKgAAAAA4QngDAAAARIo7aodFUgEAAADAEZIKAAAAIFJM1A6LpAIAAACAIzQVAAAAABwhvAEAAABao52eQRtjekp6U9IASUWSbrDWHg6z3ZmS/iApU5KVdLW1tqi5fZNUAAAAAB3DvZKWW2uHSFpeuxzOK5J+aa0dLmm8pP0t7bid9lkAAACAC9r3JWWvlXRJ7fOXJf1V0j0nbmCMGSEpzlq7VJKstaWR7JikAgAAAOgY+lpr99Q+3yupb5hthko6Yoz5izFmnTHml8aYFtsokgoAAACg/Ug1xuSfsPy8tfb5ugVjzDJJZ4R539wTF6y11hhjw2wXJ+mbks6VtFM1czBmSnqhuaJoKgAAAIBIuX+fioPW2uymVlprJzW1zhizzxiTZq3dY4xJU/i5EsWS1ltrv6x9zyJJ56uFpoLhTwAAAEDH8K6kGbXPZ0h6J8w2ayR1N8b0rl2eKKmwpR3TVAAAAAAdwxOSLjfGbJU0qXZZxphsY8wfJMlaG5D0fyQtN8ZsUE028/uWdszwJwAAACBS7g9/OmnW2kOSLgvzer6k2ScsL5V0Tmv2TVIBAAAAwJF22mcBAAAALmm/96mIGpIKAAAAAI7QVAAAAABwhOFPAAAAQKTa8UTtaCKpAAAAAOAIfRYAAAAQKZKKsEgqAAAAADhCUwEAAADAEcIbAAAAIFIMfwqLpAIAAACAI/RZAAAAQGtwR+1GSCoQVcXFxbLWul0GAAAAooimAlFjrVVMTIzKysq0Y8cOt8sBAABAlNBUIGqMMerXr5+SkpLk8/lUWlqq/fv3k1wAAID2q26itlsPj/JwaThdGGM0ePBg7d27V/v27dOOHTsUCATcLgsAAACnCE0F2kxMTIzOPvtsHTt2TJ988ok2bNigYDDodlkAAACR45KyYTH8CW2ua9euSkpKUt++fVVeXq6tW7fK7/e7XRYAAABOEn0WXNOnTx8lJyerU6dO+uSTT+Tz+ZhvAQAA0A7RVMB1/fv3V79+/fSPf/xDubm5TaQWsZLmt3Fl7UmMpAfcLsLjYiTd43YRHsbPUMtixO+h5nDhfnQg/Lg3QlMBT4iPj1fnzp01evRorVq1Sp9++mmD+RYB6UFSjCY9YqT/5Pg0a66RnuAYNelefoZaNNfwe6g5jxi3KwDgIuZUwFMSExOVmJiogQMHqqKiQoWFhQyJAgAA8DiSCnhSjx49lJSUpO7du6u4uNjtcgAAAGpw9aewSCrgaXU3zwMAAIB30WfB84xhnC4AtAfZ2dmh56mpqcrJyXGxGiBKSCrC4pAAAIBTIj8/v9n1gUBA2dnZSk9P1+LFi9uoKgBtgeFPAACgTTz77LMaPny422UAiAKaCgAAEHXFxcV67733NHv2bLdLAZypG/7k1sPFJZv4AAAgAElEQVSjaCoAAEDU3XnnnXryyScVE8OpB3A64r9sAAAQVYsXL1afPn2UlZXldinAKWFj3Xt4FU0FAACIqpUrV+rdd9/VgAEDNG3aNK1YsULTp093uywApxBNBQAAiKrHH39cxcXFKioq0htvvKGJEydq4cKFbpcF4BTy8HQPAAAAwFuskQKcQTfCIQEAAG3mkksu0SWXXOJ2GQBOMZoKAAAAIFIkFWExpwIAAACAIzQVAAAAABwhvAEAAAAiZI3kj3Xz7/JBFz+7aSQVAAAAABwhqQAAAAAiZI1RIM7NU+hqFz+7aSQVAAAAAByhqQAAAADgCMOfAAAAgFYIxMa6XYLnkFQAAAAAcISkAgAAAIiQlVFAJBUNkVQAAAAAcISmAgAAAIAjDH8CAAAAImRl5Gf4UyMkFQAAAAAcoakAAAAA4AjDnwAAAIBWCHAK3QhJBQAAAABHaLMAAACACHGfivBIKgAAAAA4QlMBAAAAwBGGPwEAAAARYvhTeCQVAAAAABwhqQAAAABagaSiMZIKAAAAAI7QVAAAAABwhOFPAAAAQISsjPwMf2qEpAIAAACAIyQVAAAAQIRqLinLKXRDJBUAAAAAHKGpAAAAAOAI2Q0AAADQCtynojGSCgAAAACOkFQAAAAAEaqZqE1S0RBJBQAAAABHaCoAAAAAOMLwJwAAACBCVuKO2mGQVAAAAABwhKQCAAAAiBh31A6HpAIAAACAIzQVAAAAABwhuwEAAAAixH0qwiOpAAAAAOAITQUAAAAARxj+BAAAALRCex3+ZIzpKelNSQMkFUm6wVp7OMx2T0r6jmoCiKWS7rDW2ub2TVIBAAAAdAz3SlpurR0iaXntcj3GmAskXSjpHEmjJI2T9K2WdkxSAQAAAESonU/UvlbSJbXPX5b0V0n3NNjGSuosqZMkIyle0r6WdkxSAQAAALQfqcaY/BMeP2rFe/taa/fUPt8rqW/DDay1qyV9JGlP7eMDa+3nLe2YpAIAAABoPw5aa7ObWmmMWSbpjDCr5p64YK21xphG8ySMMYMlDZeUUfvSUmPMN621/2iuKJoKAAAAIEJWRn4PD3+y1k5qap0xZp8xJs1au8cYkyZpf5jNpkjKtdaW1r7nfUkTJDXbVDD8CQAAAOgY3pU0o/b5DEnvhNlmp6RvGWPijDHxqpmkzfAnAAAA4FQKtN9T6CckvWWMuVXSDkk3SJIxJlvSj621syX9WdJESRtUM2k7x1r7vy3tuN0eEbQPZWVlauGyxgAAAGgD1tpDki4L83q+pNm1zwOSbm/tvmkqEFUlJSUqLy/XqlWrVF5eru3bt6tr164KBoM0GwAAAKcJmgpE1ZAhQ3To0CFNmDBBK1euVGJior7++mtVVlaGGo3NmzcrJSVFwWBQwWDQ7ZIBAACa1M7vUxE1NBVoE8YYxcTEKC0tTWlpaTp06JAuuOACrVy5Uj179tTx48dVVVWl3NxclZWVaePGjaqurtbhw4drEo3YeOkR4/bX8K6YOGkux6dZMXHSvRyjJvEz1LKYOH4PNSc23u0KALiIpgKuMsaod+/e6t27t/bu3RtqNDIyMnTw4MHQ8CkFfLrbPuJ2uZ71lHlQ823DG2LiRPPNL/SIvdvtMjzrQfMUP0MtmG9+we+hZjxlHnS7BKBNkFSER1MBzzHGqHv37urUqZNGjhypo0ePul0SACAC2dn/vB9XamqqcnJyXKwGQFuiqQAAAKdEfn5+2NcrKyt18cUXq6qqSn6/X1OnTtVDDz3UxtUBiCaaCgAAEFUJCQlasWKFkpOT5fP5dNFFF+nb3/62zj//fLdLA06Kl++o7RbuqA0AAKLKGKPk5GRJks/nk8/nkzFMegdOJyQVAAAg6gKBgLKysrRt2zbNmTNH5513ntslASelZqI2p9ANkVQAAICoi42N1fr161VcXKy8vDxt3LjR7ZIAnEI0FQAAoM10795dl156KVeGAk4zNBUAACCqDhw4oCNHjkiSKioqtHTpUg0bNszlqoCTU3efCrceXsWAMAAAEFV79uzRjBkzFAgEFAwGdcMNN+iaa65xuywApxBNBQAAiKpzzjlH69atc7sMAFFEUwEAAAC0gpeHIbmFORUAAAAAHCGpAAAAACJkZbijdhgkFQAAAAAcoakAAAAA4AjDnwAAAIAI1dynglPohkgqAAAAADhCmwUAAAC0ApeUbYykAgAAAIAjNBUAAAAAHGH4EwAAABChmonaDH9qiKQCAAAAgCMkFQAAAECESCrCI6kAAAAA4AhNBQAAAABHGP4EAAAAtIKf4U+NkFQAAAAAcISkAgAAAIhQzURtTqEbIqkAAAAA4AhNBQAAAABHyG4AAACACHGfivBIKgAAgOfMmzdPzzzzTGh57ty5evbZZ12sCEBzaCoAAIDnzJo1S6+88ookKRgM6o033tD06dNdrgqoEVCsaw+vYvgTAADwnAEDBqhXr15at26d9u3bp3PPPVe9evVyuywATaCpAAAAnjR79my99NJL2rt3r2bNmuV2OQCaQVMBAAA8acqUKZo3b558Pp9ee+01t8sBJNVM1OaO2o3RVAAAAE/q1KmTLr30UnXv3l2xsZzEAV5GUwEAADwpGAwqNzdXf/rTn9wuBUALaCoAAIDnFBYW6pprrtGUKVM0ZMgQt8sBQmruU8EpdEMcEQAA4DkjRozQl19+6XYZACJEUwEAAAC0gpfvF+EWbn4HAAAAwBGaCgAAAACOMPwJAAAAiFDNRG2GPzVEUgEAAADAEZIKAAAAoBVIKhojqQAAAADgCE0FAAAAAEcY/gQAAABEyMrIz/CnRkgqAAAAADhCUgEAAABEqOaSspxCN0RSAQAAAMARmgoAAAAAjpDdAAAAAK3AfSoaI6kAAAAA4AhJBQAAABChmonaJBUNkVQAAAAAcISmAgAAAIAjDH8CAAAAIsQdtcMjqQAAAADgCEkFAAAA0ArcUbsxkgoAAAAAjtBUAAAAAHCE7AYAAACIEPepCI+kAgAAAIAjNBUAAAAAHGH4EzzNWqtgMOh2GQAAAJIY/tQUmgp4yrFjx1RdXa3PP/9cZWVlWrVqlSorK90uCwAQgezs7NDz1NRU5eTkSJJ27dqlW265Rfv27ZMxRj/60Y90xx13uFUmgCigqYArrLU6dOiQqqqq9Nlnn6m0tFRlZWXasWOHJKlv3776+uuvdeGFF2rVqlUuVwsAiER+fn7Y1+Pi4vTUU09p7NixOn78uLKysnT55ZdrxIgRbVwhcGqQVDRGU4Goqhu+tHfvXlVVVWnt2rWqrKxUeXm59u/fr5iYGA0YMEDJycnKzc3V2WefrVWrVqlnz54yxoT2ExMfo6fMgy5+E2+LiYvRfPMLt8vwtJi4GD1onnK7DM/iZ6hlMXH8HmpOTHzT0zTT0tKUlpYmSUpJSdHw4cO1e/dumgrgNEJTgaj67LPPVFlZqdLSUsXGxmrEiBHq3LmzVq9ereHDh+vw4cPq2rVri/sJ+oKyr7dBwe2UuSko+6zbVXibuSMoS0/RJHM3P0MtMXfwe6g55qbI5r8VFRVp3bp1Ou+886JcEYC2RFOBqBo9erRyc3M1ePBg7d+/X126dHG7JACAS0pLS3X99dfrmWeeiegPSoAXWRn5Gf7UCJeUBQAAUefz+XT99dfr5ptv1nXXXed2OQBOMZoKAAAQVdZa3XrrrRo+fLjuuusut8sBHKm5pGycaw8njDHfN8ZsMsYEjTHZzWx3lTFmszFmmzHm3kj2TVMBAACiauXKlXr11Ve1YsUKjRkzRmPGjNGSJUvcLgvoiDZKuk7S35vawBgTK+k5Sd+WNELSTcaYFq+qwJwKAAAQVRdddJGstW6XAXR41trPJdW7wmYY4yVts9Z+WbvtG5KulVTY3JtoKgAAAIBWcPk+FanGmBNvCvO8tfb5U7j/dEm7TlgultTi5dpoKgAAAID246C1trn5EMsknRFm1Vxr7TvRKoqmAgAAAIhQzURt715S1lo7yeEudkvKPGE5o/a1ZjFRGwAAAECdNZKGGGMGGmM6SZom6d2W3kRTAQAAAHQAxpgpxphiSRMkvWeM+aD29X7GmCWSZK31S/qppA8kfS7pLWvtppb2zfAnAAAAIELt+Y7a1tq3Jb0d5vUSSVefsLxEUquu+0xSAQAAAMARkgoAAACgFZze2fp0RFIBAAAAwBGaCgAAAACOkN0AAAAAEfL6fSrcQlIBAAAAwBGaCgAAAACOMPwJAAAAiBDDn8IjqQAAAADgCEkFAAAA0Art9Y7a0URSAQAAAMARmgoAAAAAjjD8CQAAAIhQzURtTqEbIqkAAAAA4AhtFgAAABAhLikbHkkFAAAAAEdoKgAAAAA4wvAnAAAAoBUY/tQYSQUAAAAAR0gqAAAAgAgxUTs8kgoAAAAAjtBUAAAAAHCE4U8AAABAhKwkP8OfGiGpAAAAAOAISQUAAAAQMaMAp9CNkFQAAAAAcISmAgAAAIAjZDcAAABAhLhPRXgkFQAAAAAcIakAAAAAWoGkojGSCgAAAACO0FQAAAAAcIThTwAAAECErAx31A6DpAIAAACAIzQVAAAAABxh+BMAAAAQoZr7VHAK3RBJBQAAAABHaLMAAACAVuA+FY2RVAAAAABwhKYCAAAAgCMMfwIAAAAiVDNRm+FPDZFUAAAAAHCEpAIAAACIkJVRIEhS0RBJBQAAAABHaCoAAAAAOMLwJwAAACBSVvL7Gf7UEEkFAAAAAEdIKgAAAIAIWWsU8HMK3RBJBQAAAABHaCoAAAAAOEJ2AwAAAESoZvgTE7UbIqkAAAAA4AhJBTzP7/e7XQIAIALZ2dmh56mpqcrJyQktz5o1S4sXL1afPn20ceNGN8oDTg0rkoowaCrgWcFgUGvXrlV1dbXbpQAAIpCfn9/kupkzZ+qnP/2pbrnlljasCEBbYfgTPMdaq8LCQlVUVGjQoEFKTEx0uyQAgEMXX3yxevbs6XYZAKKEpAKeYa3VV199pbKyMg0aNEhHjhxRjx493C4LAAAgxFojv4/hTw2RVMB11lrt2bNHZWVlstYqKSlJaWlpbpcFAACACJFUwFWBQEB5eXlKSUlRUlKSBg0apL179zbaLj5WMje5UGA7ERcjmTvcrsLb4mIkc7fbVXgXP0Mti+P3ULPi+cMtOgyjYIBT6IY4InBFMBjU+vXrVVVVpXPPPVfJyclatWpVk9v7AtLv27C+9ua2oDTf7SI8bj7HqFkcn5bN5/dQs24LuF0BADcx/AltyufzafPmzaqoqFBGRoYSExOVnJzsdlkAgCi76aabNGHCBG3evFkZGRl64YUX3C4JwClEUoE2EQwGVV1drby8PPXv319JSUlKTU3Vli1b3C4NANAGXn/9dbdLAE4NK4n7VDRCUoGostbK7/crNzdX1lqdd955ysjIcLssAAAAnEI0FYgaa60KCgrk8/k0duxYJSQkKC6OcAwAAOB0wxkeosYYo+HDh2vdunXq3Lmz2+UAAAA4Zw3Dn8IgqUBUJSQkuF0CAAAAooykAgAAAIiUleQ3blfhOSQVAAAAAByhqQAAAADgCMOfAAAAgNbwu12A95BUAAAAAHCEpAIAAACIlBVJRRgkFQAAAAAcoakAAAAA4AjDnwAAAIBIMfwpLJIKAAAAAI6QVAAAAACRspJ8bhfhPSQVAAAAAByhqQAAAAA6AGPM940xm4wxQWNMdhPbZBpjPjLGFNZue0ck+2b4EwAAABApKyngdhEnbaOk6yT932a28Uu621r7qTEmRdJaY8xSa21hczumqQAAAAA6AGvt55JkjGlumz2S9tQ+P26M+VxSuiSaCgAAAOCU6SCXlDXGDJB0rqRPWtqWpgIAAABoP1KNMfknLD9vrX2+bsEYs0zSGWHeN9da+06kH2KMSZb0/yTdaa091tL2NBUAAABA+3HQWht2krUkWWsnOf0AY0y8ahqKP1pr/xLJe2gqAAAAgEid5nfUNjUTLl6Q9Lm19teRvo9LygIAAAAdgDFmijGmWNIESe8ZYz6ofb2fMWZJ7WYXSvqhpInGmPW1j6tb2jdJBQAAABCpdpxUWGvflvR2mNdLJF1d+/xjSU1fHqoJJBUAAAAAHKGpAAAAAOAIw58AAACASLXj4U/RRFIBAAAAwBGaCgAAAACOMPwJAAAAiBTDn8IiqQAAAADgCEkFAAAA0BokFY2QVAAAAABwhKYCAAAAgCMMfwIAAAAiZSX53C7Ce0gqAAAAADhCUgEAAABEykoKuF2E95BUAAAAAHCEpgIAAACAIwx/AgAAACLFHbXDIqkAAAAA4AhJBQAAABApkoqwSCoAAAAAOEJTAQAAAMARhj8BAAAAkWL4U1gkFQAAAAAcIakAAAAAWoOkohGSCgAAAACO0FQAAAAAcIThTwAAAECkmKgdFkkFAAAAAEdoKgAAAAA4wvAnAAAAIFIMfwqLpAIAAACAIyQVAAAAQKSsJJ/bRXgPSQUAAAAAR2gqAAAAADhCUwEAAE6J7Ozs0OOqq66qty4nJ0dnnXWWBg8erCeeeMKlCoFTwEoKuPjwKOZUwPOqq6vdLgEAEIH8/PywrwcCAc2ZM0dLly5VRkaGxo0bp8mTJ2vEiBFtXCGAaKGpgKdt3bpVfj/XbQOA9iwvL0+DBw/WoEGDJEnTpk3TO++8Q1OB9otTk0YY/gRPstaqsrJS1dXVSkxMdLscAIADu3fvVmZmZmg5IyNDu3fvdrEiAKcaTQU8acOGDZLEX7EAAADaAYY/wVMCgYDKy8t1xhlnqLS0VMYYt0sCADiUnp6uXbt2hZaLi4uVnp7uYkWAA9xROyySCniGtVaffvqp4uLiQuNuAQDt37hx47R161Z99dVXqq6u1htvvKHJkye7XRaAU4ikAp5QXV2t8vJyDRw4UEVFRW6XAwA4CdnZ2aHnqampysnJkSTFxcVpwYIFuvLKKxUIBDRr1iyNHDnSrTIBZ0gqwqKpgOsqKyu1bt06JSQkqF+/fjQVANBOrVmzpslhq1dffbWuvvrqNq4IQFuhqYCrgsGgPv30Uw0bNkxffPFFo/XWWsXExCi5SxfdVlHhQoXtQ4yk+W4X4XEco+ZxfFoWK+k2t4vwsDhj1L9/f/Xp06fRuhNTCwCnJ5oKuKa0tFTl5eU677zz1K1bt7Db+P1+jR07ViX79jV6b2FhoUaNGsUlZ09grdW6des0dOhQJScnu12OpxQVFalLly7q27ev26V4yuHDh3XgwAENHTrU7VI8pbKyUhs2bFBWVpZiYph+WKeqqkqfffaZhgwZou7du9db5/f7NWfOHGVmZurxxx/nuOH0ZSX53C7Ce2gq4IpAIKCCggIlJiaGbSjqEoq8vLyw762srFSXLl20adOmtii33fD5fAoEAvr888/dLsVzKioqlJCQoJ07d7pdiqdYa1VRUaGjR4+6XYrnVFVVafXq1erUqZPbpXhKMBjU+vXr1blzZ8XGxoZev/fee3X06FGtXr1af/zjH3XGGWc0ei+JBXD6oqlAmzt06JAqKip00UUXaf369WG3qUsoGiorK9OmTZs0btw4JSUlRbvUdsXv9+vTTz/VuHHjFB8f73Y5npOXl6dx48ZxmeIw8vLylJ2dzbFpIBAIaO3atTrnnHOUkJDgdjmeUpdYDB06NPSHoeXLl0uqOW4/+9nPFB8frwULFpBY4PRjJQXcLsJ7+C8dbWr//v3asmWLkpKSmhy2ZK2VtbbR63UNxciRI2kowigqKlJGRgYNRRg+n09xcXGcNDchMTFR5eXlbpfhObGxsRo0aJC2b9/udimek5CQoLPPPlubN2/WsWPH6q2LjY3Vr371Kxlj9KMf/UiBAGdfQEdAUoE24/P59NVXXyk7O1tr1qxpcpvY2NhG64PBoCoqKtS5c2cVFha2RbntSjAYVGVlpRITE7V79263y/Ecv98vv9/f5M9dR1ddXa2CggIa0iaUl5crNze33lAf1Ki7v9CJQ6HqhkFJNX8M6tu3r/r371+vqWcYFHD6oalAm9i5c6d8Pp8uvPBCxcU1/rGrSyestY2GPZWXl2vjxo3KyspSSkpKW5XcblhrVVBQoGHDhjU54b2j27FjhxISEsKO8Yb09ddf69ChQxoyZIjbpXhSeXm5CgsLlZWVRdoVRkVFhTZs2KBhw4YpJSUlNAxKqvmDx2OPPaadO3fqlVdeoXHF6YP7VDRCU4GostaqqqpKhw4dUmJiYpMNhd/vbzGhCHfJWfzzr/BbtmxxuxTPqpukvWvXLrdL8aS6ydpHjhxxuxTPqqqq0qpVq5i03QRrrdauXdvk5G2fz6czzzxT/fr1a9SYkVoApweaCkSNtVZbt25VMBjU6NGjlZubG3Y7v98fNqGo++vX2LFj1bVr17Youd2pm0ialZXFyU4zmKTdMiZrN6/uQghjxozhr+1NqEuVhw8fHrqkdV1qYa3Vb37zG61Zs0ZvvfUWE9/RvnFH7bBoKhBVKSkp6tKlS9irfzR32dgTE4rNmze3RantUlVVlYwxKigocLsUz6r7K3x+fr7bpXhaZWWl8vLyuFJPM/x+v3Jzc9W5c2e3S/Esa63y8/Mb/d6vSyz8fr/69++vfv36NfpZI7EA2jeaCkSNMUZpaWnasWNH2PVNXTa27lKF5557LnMEmlFRURGaa8KJYNO4uVtkuDlgy6y1Wr9+vQYPHsz8rmbUXalvxIgRoSv1nZhYvPDCC1qyZIkWLVrElfyA0whnImhzJ07Kbijctc8R3rZt2zR48GAaihYcP36cE8AIpKSk6Pjx426X4WnGGA0ZMkRbt24N+/uro/j666917bXXasyYMbr22mt1+PDheuuTkpI0YsQIPf300xo9erTGjBmjP/7xj5JqjuHYsWO1adMm9e7dWz/+8Y9Dx/JPf/qTRo4cqZiYmEbJ4uOPP67BgwfrrLPO0gcffBB6PScnR2eddZYGDx6sJ554IsrfHKhVd0dttx4eZVr4xdhxf2vCsUAgoOrqauXm5uqCCy7QqlWrdP7554fGJje8dvmJQ564dGPz/H6/fL7/397dB1VZp38c/3w5BIoICGKArhAPamrlA7SRJT0QPlRkzf5MS1sXG91qWlwjV2uz2d2wmmp7sjFbS6lfqxm16aKb6U5OjRm4RlBhKyISoJhCgojAOZz790fJL7dTHvdWD+j7NcMfHs65+XZPMV1+ruu+nOrZs6evj9LlHT16VAEBAfw7dQLH2sR+bH8M/l9ra6scDsc5O1uxdOlShYSEaOrUqVq5cqUOHz6sWbNmHfeepqYm3XXXXXrqqafUs2dP3XPPPVqyZIlyc3NVWlqqfv36yeFwqLq6WrGxserTp49aW1slffu0tuHDh2vr1q2SpLKyMk2dOlVFRUXau3ev0tPTOx9MMWjQIG3cuFEDBgxQSkqKVq5cqaFDh57ZG4JTqVsMdZnIZEuTfNhSu8xstywr2XcH8Iz2J5xRLte3k03/2fZ07Dn5l1xyifr06eOLo3Ubbrdb27dv18iRI+nt9sKxAWQSnRNjWNs7TqdTxcXFGjlypMcn2p3tZs+erfXr1ysgIEBlZWXatm2bvv76a61YsaLz9/ebb76pCRMm6KqrrtKzzz6r5uZmzZo1S3PnztWhQ4e0fft2FRcXa9KkSaqqqlJ6erpefPFFGWM0ZswYuVwuJSUlKS4uTqmpqZoyZYqee+45vf766zpw4ICSkpJUVVWltLQ0xcfHKy4uTq2trUpPT1dMTAwzVIAPnHu/DeETx4ayPS0fsyxLLS0tCgwM1K5du3xwuu6lvb1dlmXps88+8/VRurxjf/u+fft2Xx+lW2BY23sul0tbt249Jwv7ffv2qbq6WkuXLtWFF16okpISJSQkaN68eZ2JxbEHcBQVFenVV19VRkaGgoOD9cQTT6h///7atm2bbrzxRgUHByswMFDLly/Xxo0bFR4ertLSUoWFhWngwIGqqqpSfX295s6dq2nTpun+++/XzJkzFRkZqYKCAsXHx3eea+HChSorK9PixYt9cl9wDrEksSj+BygqcEY4nU6PQ9lOp1MlJSW6+OKLFR4e7oOTdS9tbW0qKSnR6NGjaefxwqFDh7R//34NHjzY10fpFiorK9WrVy/169fP10fp8o5tkh4yZMhZOWycmZmp/fv3/+D1hQsXyt/fXykpKZ2JxVtvvaVZs2YpLS1N27ZtU2xsrC6//HIZY9TQ0KAJEyYoIiJCTU1N6t+/v3bs2KGsrCz17t1beXl5evrppxUWFqb8/HxZliW32628vDxNmDBBZWVluuSSS5STk6O8vDytXr1akrRu3TqVl5dr3759+te//nXcjpV3331X2dnZ6ujo0J133qn58+efsfsGnMsoKnBaOZ1Oud3uH8xPSN+28ZSUlCg2NlahoaEe34PjVVRUKC4uTpK4X15obGxUr169uFde6tWrlxobGxUREeHro3QL8fHxKi8v10UXXeTro5xyf/vb3370e5GRkaqtrdXXX3+tjo4O9e3bV6+99lpna+bTTz+twsJCRUVFKSAgQDExMaqpqdGQIUP0+uuvKzg4WM8995xuvfVWffnll4qKitLWrVuVmJioXr16ac+ePdq5c6cmTJig5cuXS5Kys7NljNFjjz2mqqoqVVRUaMSIEQoLC9OGDRt0wQUXKDc3V5Zl6Y033lBhYWHnnEVmZiZzFji12FPhEYPaOG2ObdMuKir60f+pc7vdtFqcBO7XyTn2+40ZAe9wv07e2f7fZE5OjhoaGo57raGhQT//+c/10UcfaerUqTp8+LC2bNmihoYGFRQUqL6+XtnZ2XK73Ro/frycTqf++c9/6rbbblNBQYG++eYbzZgxQ8uXL1doaKjcbreqq6sVGBio6Oho7dq1S2FhYUpISNDnn38up9OpwMBApaSkqLq6Ws3NzUpISFB4eLh27NjR+bjaSZMmafbs2frd736ndevWaezYsbrZ9OUAABNXSURBVHr00UclSQsWLDjj9w7/lW7xy8dEJFu63odzO68xqI1zjDFGPXr00NixY319FADAf8HTPFJ9fb0mT56s9vZ2ffnll1qzZo0uuOAC9e3bV6+88orq6upUV1enIUOGqKCgQO3t7XrmmWdUX1+vkJAQ3XHHHXr++ed18OBBNTc3a968eXrrrbcUExOjTz/9VNK3rYsHDhxQQECAYmNj9Ytf/EJ//vOf1draqtDQUF1++eVasWKFHA6HEhIS5HA4FB0drZaWFiUlJamoqEh33HGHOjo6dPToUb311lsMbwOn2dn71ysAAOCUSU9P1/Dhw5WWlqb9+/crJCRE5eXl+vDDD9Xe3q5bb71V0dHRuvbaaxUSEqLp06dr6tSpCg8P16RJk9TS0qJdu3ZpypQpCgoKUnR0tKZNm6Zly5Zp7969mjNnjuLi4pSenq5BgwZp6NCham1t1U033SSn06mHH35YoaGham9vV2trq6qqqvTOO+9oxYoVGjhwoJYuXaolS5aotrZWw4cPl/RtQjFlyhQKCpx6Lh9+dVEkFQAA4IQ2bdp03J+PJRY5OTny8/PTjBkzlJGRIbfbraamJjU1NWnDhg166KGHlJKSosOHD+uiiy7StGnTVF9fL4fDoYKCAjU3N8vtduvDDz/Up59+qvDwcDmdTjU2NsqyLCUkJOjee+9VVFSUOjo6dN1112nTpk06cOCAfvvb36qurk7V1dWaM2eOqqqqFBcXp+XLl6u2tla5ubm6/fbbfXTHgHMLMxUAAMCW+++/X3369NFDDz2k++67Tx0dHdq0aZN2797duam9oaFBSUlJio+PV0ZGhl566SVdccUVGjVqlF544QWdd955SkpKUkdHh3bu3KmVK1cqKytLjY2Ncjqdmjdvnp577jm9//77uvzyyzVu3Djdcsstio+P1zXXXKOMjAyFhYVp7dq1ysjI0EcffaS9e/eqZ8+eeuKJJ36woA9dUveYqQhPtnStD9Ov/K45U0FRAQAAvJaenq66urrjXuvo6FB7e7t2796ta6+9VqtXr9Z9992nFStWqHfv3nrggQc0f/58DR8+XNXV1fL399cjjzyiOXPmKDg4WO3t7fLz85PT6VRWVpY2bNigxsZGHT58WP7+/mppadHVV1+tfv36aePGjaqvr5ckDR8+XEFBQWpsbNS///1vDR48WAcOHFCPHj3kcDh0zTXXKC8vT4GBgRowYICmT5+uhQsX+uK2wTsUFd7ookUF7U8AAMBr/9kGdUx+fr5mzpyp1157TaGhofr73/+u3r176+DBg52PdpWk3NxcffbZZ7IsS4GBgbrhhhtUXFzcWTzcfffd+utf/6rbbrtNn376qUpLS5WQkKAVK1Zo9OjRSktLU3FxsWpqarRq1Sq5XC6lp6crISFBZWVlCgsL09y5c/XAAw9o8+bNWr16te655x7169dPOTk5Z/JWAecUBrUBAMApERsbq7y8PBUVFcnf31+JiYkKCAjQlClTtGbNGklSWlqa3nvvPeXn56tXr156//33lZycrPLycvXu3VvPP/+84uPjtWXLFo0dO1YhISFKSUlRbm6u+vXrp9raWk2YMEHBwcHKy8vTLbfcoqVLlyowMFDGGPXs2VPV1dWdbVeWZen999/vHN4GbDu2UdtXX10U7U8AAMC2rVu36sEHH5QxRmVlZTpy5Iiys7P1pz/9SdOnT9fq1avlcDj0y1/+UikpKfr1r38tl8ulZcuWqba2Vo888oj69u2r/fv3yxijm266SXPnztWECRPUo0cPHTp0SJZlaeDAgZ3L7dxut1wul/r06aO2tjYNGjRIpaWlCg4OVmtrq9xut9rb23X++edrxIgRioiIUGFhoYKCgrRixQqNGjXK17cNx+se7U99ki1d7cP2p791zfYnkgoAAGBbSkqK9uzZo2XLlunpp5+Wv7+/pkyZoo6ODr333nuaPHmyHn30Ua1evVqXXXaZzj//fE2cOFFZWVlqb2+Xy+VSUVGRUlNT5XA49MYbbyg4OFjNzc1asGCBnnjiCXV0dOjtt99WWFiYLMvSpZdeqqNHj8qyLM2ZM0ejRo1SRESEUlJSdPToUS1ZskQ333yz6urqlJqaqnfeeUfl5eV66aWXdNddd/n6lqG7OrZRm0fKHoeZCgAAYJu/v78WL16scePG6ciRI4qIiNCwYcOUlZUlh8OhoUOHavbs2Vq2bJmuvPJKuVwuzZgxQ5JUWFio8847T1deeaWampokSX5+flq3bp2ioqL01FNPye12KzAwUEeOHFFRUZGio6NVWVmpoUOHKiYmRlu2bFFISIjS09NVUVGhxMREhYeHa9WqVZKkvXv3KigoSPX19brssst06NAh7du3T9HR0b66ZcBZhaQCAACcEhMnTtTOnTtVVVWljo4OVVZWKiMjQ62trcrMzNTmzZt18OBBuVwuXXXVVZ1D3zt27FB6erp27Nih5ORktbW1adSoUdq2bZucTqf27Nmju+++W62trZo0aZLWrl2rxsZGbdiwQYsXL9ZXX32lDz74QKNHj5a/v78aGhq0b98+VVRU6Oabb1ZMTIzy8/Pl5+enzz77TKGhodq7d6/S0tL0xz/+0cd3DTg7MFMBAABOufXr12vOnDlqbGzUz372MxUWFioyMlKTJk1SYGCgtmzZoujoaFVUVGjv3r3auHGjSkpKVFpa2jl/0djYqOzsbP3qV7/S1KlTFRkZqZKSEjU0NGjEiBEqKirSoEGDtHv3bvXu3VtOp1P9+/fX9OnTNX/+fI0dO1aVlZVqa2uTy+XSk08+qQsvvFBPPvmkjh49qscff1zJyV2uNf1c1j1mKkKTLY3x4UzFP5ipAAAA54hjqcU777yjiIgIFRUVKSUlRUlJSRo4cKBGjx6tTz75RJI0YMAA+fn5ac2aNbr99tvl5+enyspKud1uvfLKKxo3bpzGjRsnh8OhtWvXKiYmRl988YXi4+MVFhamRYsWacGCBRo2bJjq6ur05ptvKjAwUI8//rhGjhwpSQoNDdUf/vAH5eTkyOl0qqamRjExMfrNb36jxMREXXzxxZ3nAXDyKCoAAMBpk5KSovLycn3yySfq37+/Vq1apeuvv17r16/X+PHjVVZWppaWFj3zzDOqra1VaWmprrnmGr366qvq0aOHoqKiNGTIEG3cuFHl5eUKDg7WkSNHFBISoqysLJWXl2vlypXKzMxUZGSk/Pz8VFFRobi4OE2bNk3nnXeeoqKilJGRoZEjR2rkyJHavHmzqqurlZmZqeLiYoa3cXIsSU4ffnVRFBUAAOC0OTbA/cgjjyg/P1+TJ0/WkSNHZFmWGhoaFBAQoNmzZ6usrEw7d+7UX/7yFz322GNas2aNQkNDddNNN6mwsFAlJSVavHixCgoKFB8frzvvvFMvvviiWlpalJqaqmHDhum2225TUFCQLMvSgQMHlJiYqPb2djU3N2vRokWKj4/Xu+++K0n64IMPFBkZqfLychljjhveBnDyKCoAAMBpNXHiRL399ttKTU3Vgw8+qNraWsXGxuqKK66QJMXFxSktLU1XX321lixZovj4eNXW1qqlpUW5ublasGCB+vTpo0svvVS1tbVqa2vTjTfeqLfffltRUVG67rrrJEk1NTW64447FBQUpJycHOXn52v79u06cuSIYmJi9MILL+irr75Sz549lZycLIfDIUk6ePCgpG/bsGpra31zk4BujqICAACcdsfaoCorK+VyuVRRUaHMzMzj3pOZmam8vDxJUmNjo8aMGSNjjDIzM9Xc3Ky2tjY1NTVp//79uvTSS5WSkqLGxkbt379f7e3tWrVqldxut5xOpx5++GHl5+frhhtuOO5nHFuuJ6lzoV5ERMSZuQk4e7BR+wcoKgAAwGn3/T0W9913X+cei4ULF2rjxo3q37+/Zs6cqfr6eiUmJqqpqalzj8XgwYPl7++vK6+8Uv/4xz80ceJEORwO+fv7a9CgQXrsscd04YUXasyYMVq7dm3n4PeqVas0depUnX/++Z1tTS+//LJaWlp0ySWXqLy8XHfffXdnkVFTU6P+/fv76hYB3RpFBQAAOCM87bH4/e9/r5KSEmVmZqpHjx568803tWvXLi1atKhzj0V+fr6uv/567d69W5s3b1ZxcbHa2tpUWVmpb775Rnv27FFFRYWOHj2qQ4cOac+ePRoxYoQSExM1fvz44xIQh8Oh7OxslZSU6JVXXtHHH38sy7L08ccfKzQ0lGV4OLFuvFHbGPM/xpgvjDFuY8xPPpbWGOMwxhQbYwq8uTYbtQEAwBn1/dSio6NDWVlZnalFcnKyMjMzNXPmTE2fPv0Hm7GHDRumyZMna+jQofL399cLL7zQORvR2toqh8Mhl8ulgwcPKjU1VZI0f/58TZ48WS+//LJiY2O1evVqSd8WOevXr1diYqKCgoK0fPly39wQ4Mz5XNItkpZ68d5sSTskhXhzYZbfAQAAoCvoHsvveidbGunD5Xcf2l9+Z4zZLCnHsiyP/yDGmAGS8iTlSpprWdYNnt73fSQVAAAAgLeOtT/5Tl9jzPeLgZcsy3rpFP+MZyTNk9Tb2w9QVAAAAADdx8GfSiqMMZskRXn41oOWZa050cWNMTdI+tqyrO3GmKu8PRRFBQAAAHCWsCwr3eYlxkjKNMZMlNRDUogx5n8ty5r2Ux+iqAAAAAC8ZUly+voQp49lWQskLZCk75KKnBMVFBKPlAUAAADOCcaYm40xNZJSJa0zxmz47vUYY8x6W9fm6U8AAADoArrH05+Cki0N8eHTn4rtP/3pdCCpAAAAAGALRQUAAAAAWxjUBgAAAE6Gb/dUdEkkFQAAAABsIakAAAAAvOX7jdpdEkkFAAAAAFsoKgAAAADYQvsTAAAA4K2zfKP2f4ukAgAAAIAtJBUAAACAtyxJHb4+RNdDUgEAAADAFooKAAAAALbQ/gQAAAB4iz0VHpFUAAAAALCFpAIAAAA4GSQVP0BSAQAAAMAWigoAAAAAttD+BAAAAHiLjdoekVQAAAAAsIWiAgAAAIAttD8BAAAA3rIkdfj6EF0PSQUAAAAAW0gqAAAAAG+xUdsjkgoAAAAAtlBUAAAAALCF9icAAADAW7Q/eURSAQAAAMAWkgoAAADAW2zU9oikAgAAAIAtFBUAAAAAbKH9CQAAADgZbNT+AZIKAAAAALaQVAAAAAAnw/L1AboekgoAAAAAtlBUAAAAALCFogIAAACALRQVAAAAAGyhqAAAAABgC0UFAAAAAFsoKgAAAADYQlEBAAAAwBaKCgAAAAC2sFEbAAAA8JolyenrQ3Q5JBUAAAAAbKGoAAAAAGAL7U8AAACA1yxJLl8fosshqQAAAABgC0UFAAAAAFtofwIAAAC8xtOfPCGpAAAAAGALSQUAAADgNQa1PSGpAAAAAGALRQUAAAAAW2h/AgAAALzGoLYnJBUAAAAAbCGpAAAAALxGUuEJSQUAAAAAWygqAAAAANhC+xMAAABwUthT8Z9IKgAAAADYQlIBAAAAeI1BbU9IKgAAAADYQlEBAAAAwBbanwAAAACvWWJQ+4dIKgAAAADYQlIBAAAAeI1BbU9IKgAAAADYQlEBAAAAwBbanwAAAACvMajtCUkFAAAAAFtIKgAAAACvMajtCUkFAAAAAFsoKgAAAADYQvsTAAAA4DUGtT0hqQAAAABgC0UFAAAAAFtofwIAAAC8xtOfPCGpAAAAAM4Bxpj/McZ8YYxxG2OSf+J9YcaYfGPMl8aYHcaY1BNdm6QCAAAA8Fq3HtT+XNItkpae4H3PSnrXsqxfGGMCJAWd6MIUFQAAAMA5wLKsHZJkjPnR9xhjQiWNlTTju8+0S2o/0bVpfwIAAAC6j77GmH9972vWKb7+BZIOSFpujCk2xiwzxvQ60YdIKgAAAACv+XxQ+6BlWT81D7FJUpSHbz1oWdYaL67vL2mUpHstyyo0xjwrab6kh070IQAAAABnAcuy0m1eokZSjWVZhd/9OV/fFhU/iaICAAAAOCnddlD7hCzLqjPGVBtjBluW9W9J10oqO9HnmKkAAAAAzgHGmJuNMTWSUiWtM8Zs+O71GGPM+u+99V5JrxtjSiWNkLTohNe2LOunvv+T3wQAAABOkR9/JFEXYkySJf3ZhyfI3P5TMxW+QvsTAAAA4DWfD2p3SbQ/AQAAALCFpAIAAADwGkmFJyQVAAAAAGyhqAAAAABgC+1PAAAAgNcsnc17Kv5bJBUAAAAAbCGpAAAAALzGoLYnJBUAAAAAbKGoAAAAAGAL7U8AAACA1xjU9oSkAgAAAIAtFBUAAAAAbKH9CQAAAPAaT3/yhKQCAAAAgC0kFQAAAIDXGNT2hKQCAAAAgC0UFQAAAABsof0JAAAA8BqD2p6QVAAAAACwhaQCAAAA8BqD2p6QVAAAAACwhaICAAAAgC20PwEAAABeY1DbE5IKAAAAALaQVAAAAABeY1DbE5IKAAAAALYYy7J8fQYAAACgWzDGvCuprw+PcNCyrPE+/PkeUVQAAAAAsIX2JwAAAAC2UFQAAAAAsIWiAgAAAIAtFBUAAAAAbKGoAAAAAGALRQUAAAAAWygqAAAAANhCUQEAAADAFooKAAAAALb8HyFAlTFZw+xZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x864 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pp.plot_grid(g, cell_value=u[1::2], figsize=(15, 12))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To understand the inner workings of the discretization, and to recover the traction on the grid faces, some more details are needed. The MPSA discretization creates two sparse matrices \"stress\" and \"bound_stress\". They give define the discretization of the cell-face traction:\n",
    "\\begin{equation}\n",
    "T = \\text{stress} \\cdot u + \\text{bound_stress} \\cdot u_b\n",
    "\\end{equation}\n",
    "Here $u$ is a vector of cell center displacement and has length g.dim $*$ g.num_cells. The vector $u_b$ is the boundary condition values. It is the displacement for Dirichlet boundaries and traction for Neumann boundaries and has length g.dim $*$ g.num_faces.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The global linear system is now formed by momentuum balance on all cells. A row in the discretized system reads\n",
    "\\begin{equation}\n",
    "-\\int_{\\Omega_k} f dv = \\int_{\\partial\\Omega_k} T(n)dA = [div \\cdot \\text{stress} \\cdot u + div\\cdot\\text{bound_stress}\\cdot u_b]_k,\n",
    "\\end{equation}\n",
    "\n",
    "The call to mpsa_class.assemble_matrix_rhs(), creates the left hand side matrix $ \\text{div} \\cdot \\text{stress} $, the right hand side vector, which consists of $\\text{f}$ and $-\\text{div} \\cdot \\text{bound_stress}$ (note sign change).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also retrieve the traction on the faces, by first accessing the discretization matrices stress and bound_stress"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1MAAAKaCAYAAADbKANUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xl0VPXh///XzUwgZAj7FpIAUpVVQEAQDSoKqK2iKCoqqEVsP596TqHWr61L3QpaeypqpfZz2qpVqFKt1SpKBIWqCCGEVQRxY8lCgbBmzyz39we/pAKZhFzIvXfufT7O6cEwk+GVd+fe3Ne83/MewzRNAQAAAACaJsnpAAAAAACQiChTAAAAAGABZQoAAAAALKBMAQAAAIAFlCkAAAAAsCDYyO1s9QcAAADAz4x4NzAzBQAAAAAWUKYAAAAAwALKFAAAAABYQJkCAAAAAAsoUwAAAABgAWUKAAAAACygTAEAAACABZQpAAAAALCAMgUAAAAAFlCmAAAAAMACyhQAAAAAWECZAgAAAAALKFMAAAAAYAFlCgAAAAAsoEwBAAAAgAWUKQAAAACwgDIFAAAAABZQpgAAAADAAsoUAAAAAFhAmQIAAAAACyhTAAAAAGABZQoAAAAALKBMAQAAAIAFlCkAAAAAsIAyBQAAAAAWUKYAAAAAwALKFAAAAABYQJkCAAAAAAsoUwAAAABgAWUKAAAAACygTAEAAACABZQpAAAAALCAMgUAAAAAFlCmAAAAAMACyhQAAAAAWECZAgAAAAALKFMA4COtW7fWt99+W+9tf/3rX5WdnW1zIgAAEhdlCgBOsVdeeUXDhw9X69atlZ6erssvv1zLly+3/HiGYejrr78+6u9KS0t11113qVevXgqFQurRo4cmTZqkVatWNfhYZWVl6t27t6UcNTU1evjhh3XGGWcoFAqpV69emjZtmrZv327p8ZrLww8/rClTpjgdo1lt375dhmEoEok4HQUAfI0yBQCn0Jw5czRz5kzdd9992r17t3bu3Kmf/OQn+te//tXkx4p3oVxdXa2LL75Yn332mRYuXKjDhw9ry5Ytmjx5shYtWtSkx2qKSZMm6e2339Yrr7yiQ4cOacOGDRo2bJg+/PDDk35snHoULQCwgWmaDf0PAHCCDh48aIZCIfO1116Le59Vq1aZ5557rtm2bVuzW7du5p133mlWV1fX3S7JnDt3rnn66aebvXr1MkePHm1KMlNTU81QKGQuWLDA/POf/2x269bNLCsrazDPsY9V+3dfffWVaZqmWVJSYl555ZVmWlqaec4555gPPPCAef7559f7WEuWLDFTUlLMnTt3xv33ioqKzCuvvNJs3769+b3vfc/805/+VHfbQw89ZE6aNMm8+eabzdatW5sDBw40t27daj722GNm586dzczMTPP999+vu/+FF15o/vKXvzTPOeccMy0tzZwwYYK5b98+0zRNc9myZWZGRsZR/3bPnj3NJUuWmIsWLTKTk5PNYDBohkIhc9CgQaZpHvn/Ztq0aWa3bt3M7t27m/fff78ZiUTq/TkqKirMW265xWzXrp3Zt29f84knnjjq3ysqKjKvueYas1OnTmavXr3MZ555pu62qqoqc8aMGWZ6erqZnp5uzpgxw6yqqjoq9xNPPGF27tzZ7Natm/nmm2+a7777rnnGGWeY7du3N2fPnl33WNFo1Hz88cfN3r17mx06dDCvu+66ujHIysoyJZmhUMgMhULmihUrzBdffNE877zzzJkzZ5odOnQw7733XrN9+/bmxo0b6x5z9+7dZqtWrcw9e/bE/f8RAHCcuH2JmSkAOEVWrlypqqoqTZw4Me59AoGAnnrqKZWUlGjlypX68MMP9dxzzx11n7feekurVq3S5s2b9fHHH0uSNmzYoLKyMt1www364IMPdOmllyoUCjWa6buPdaw777xTKSkp2rVrl1544QW98MILcR/ngw8+0IgRI5SVlRX3PpMnT1ZmZqaKi4v1j3/8Q/fdd5+WLl1ad/s777yjqVOn6sCBAzr77LN16aWXKhaLqaioSA8++KB+/OMfH/V4L7/8sl544QXt2rVLwWBQP/3pTxv9eS+77DLdd999uuGGG1RWVqYNGzZIkm677TYFg0F9/fXXWrdunRYvXqy//OUv9T7GI488ou3bt+vbb7/VkiVLNH/+/LrbYrGYrrzySg0ePFhFRUX68MMP9fTTT+v999+XJM2ePVu5ublav369NmzYoLy8PM2aNavu+//zn/+oqqpKRUVFevTRR3XHHXdo/vz5WrNmjT755BP9+te/1rZt2yRJzz77rN566y199NFHKi4uVvv27XXnnXdKUt3z4uDBgyorK9OoUaMkSatWrVLv3r21e/du/epXv9LkyZOPyv/qq6/qkksuUefOnRsdSwDACWioaTnS+wAgQc2fP9/s2rVrk77nqaeeMq+++uq6ryWZH3744VH30Xdmk0zTNC+55BLzF7/4Rd3X69atM9u2bWumpaWZZ5555gk9ViQSMYPBoLlly5a62+699964M1PTp083b7jhhrg/x86dO82kpCTz8OHDdX/3y1/+0rz11ltN0zwyMzV27Ni6295++20zFArVzQ4dPnzYlGQeOHDANM0jM1Pf/Rk///xzMzk52YxEIg3OTNX+WzfffHPdbf/5z3/MFi1amBUVFXV/98orr5gXXXRRvT/LaaedZubk5NR9/ec//7nu38vNzTWzsrKOuv9jjz1m3nbbbaZpmmbv3r3Nd999t+62nJwcs2fPnqZpHpmZSklJOe5nzs3Nrbv/0KFDzTfffNM0TdPs27ev+cEHH9TdVlxcbAaDQTMcDpvbtm0zJZnhcLju9hdffPG4bLV5Y7GYaZqmOWzYMPPvf/97vT83ACCuuH0p6GCPAwBP6dixo0pKShSJRBQM1n96/fLLL3XXXXcpPz9fFRUVikQiGjZs2FH3aWj2p/bf2bVrV93XQ4YM0cGDB/XBBx9o+vTpJ/RYe/fuVSQSOer2nj17Nvhvfvnll3FvLy4uVocOHZSWlnbU4+Xn59d93bVr17r/btWqlTp16qRAIFD3tXRkg4x27dodl71nz54Kh8MqKSmJmyGeHTt2KBwOKz09ve7vYrFY3LEpLi4+6rbv/veOHTtUXFxcl1GSotGoRo8eXfe93x3Hnj17qri4uO7rjh07HvczHzsuZWVldf/WxIkTlZT030UkgUBAu3fvjvuzHvszjRw5Uqmpqfr3v/+t9PR0ff3115owYULc7wcANA3L/ADgFBk1apRatmypt956K+59/vd//1d9+/bVV199pcOHD+uxxx6TaZpH3ccwjAb/nUsuuUSLFy9WeXl5o5niPVbnzp0VDAZVUFBQ93c7d+6M+zhjx45VXl6eCgsL6729e/fu2r9/v0pLS496vIyMjEYzxnNstuTkZHXq1EmhUEgVFRV1t0WjUe3du7fu62N/5qysLLVs2VIlJSU6ePCgDh48qMOHD+vzzz+v999NT08/6uf8bo6srCyddtppdY9z8OBBlZaW6r333pN0ZBx27NhxVO7u3btb+vmzsrK0aNGio/6tqqoqZWRkxP3/tb6/v/XWWzV//nzNmzdPkyZNUkpKiqU8AIDjUaYA4BRp27atHn30Ud1555166623VFFRoXA4rEWLFumee+6RdGRL8zZt2qh169b64osv9Mc//rHRx+3atetRnw11yy23KD09XRMnTtSmTZsUjUZVVVV11CxQYwKBgK655ho9/PDDqqio0ObNm/XSSy/Fvf/YsWM1btw4TZw4UWvWrFEkElFpaan+7//+Ty+88IKysrJ03nnn6d5771VVVZU2btyo559//qS2KJ8/f742b96siooKPfjgg5o0aZICgYDOPPNMVVVV6d1331U4HNasWbNUXV1d931du3bV9u3bFYvFJB0pR+PHj9fPf/5zHT58WLFYTN98840++uijev/d66+/Xo8//rgOHDigoqIizZ07t+62ESNGKC0tTU888YQqKysVjUa1adMmrV69WpJ04403atasWdq7d69KSkr06KOPWh6D//mf/9H9999fV8727t1btytk586dlZSUFPczw75rypQpevPNNzV//nzdcsstlrIAAOpHmQKAU+jnP/+55syZo1mzZqlz587KysrS3LlzdfXVV0uSfve73+mVV15RWlqa7rjjDt1www2NPubDDz+sW2+9Ve3atdNrr72mlJQULVu2TP3799cPfvADtWnTRn369NHq1av12muvnXDWuXPnqqysTN26ddNtt92mH/7whw3e/x//+Ie+//3v64YbblDbtm01cOBA5efna+zYsZKObG6wfft2de/eXRMnTtQjjzxSd5sVU6dO1W233aZu3bqpqqpKv//97yUdKa3PPfecpk+froyMDIVCIWVmZtZ933XXXSfpyJK6oUOHSjqymUVNTY369++v9u3ba9KkSUctlfyuBx98UJmZmTrttNM0duxYTZo0SS1btpR0pIQuXLhQ69ev12mnnaZOnTpp+vTpOnTokCTpgQce0PDhwzVo0CCdddZZGjp0qB544AFLP/+MGTM0YcIEjR8/XmlpaTr33HPrPkcsNTVV999/v84//3y1a9dOubm5cR8nKytLQ4cOlWEYdcsRAQCnhnHs8pJjNHgjAADN4aKLLtKUKVOOew+YE/74xz9qwYIFcWeyEsG0adPUvXv3o3YWBACcsLjr75mZAgDgO3bt2qVPP/1UsVhMW7du1ZNPPtngdvdut337dv3zn//U7bff7nQUAPAcyhQAAN9RU1OjH//4x0pLS9PFF1+sq666Sj/5yU+cjmXJr371Kw0cOFD/7//9P5122mlOxwEAz2GZHwAAAADExzI/AAAAADiVKFMAAAAAYAFlCgAAAAAsoEwBAAAAgAWUKQAAAACwgDIFAAAAABZQpgAAAADAAsoUAAAAAFhAmQIAAAAACyhTAAAAAGABZQoAAAAALKBMAQAAAIAFlCkAAAAAsIAyBQAAAAAWUKYAAAAAwALKFAAAAABYQJkCAAAAAAsoUwAAAABgAWUKAAAAACygTAEAAACABZQpAAAAALCAMgUAAAAAFlCmAAAAAMACyhQAAAAAWECZAgAAAAALKFMAAAAAYAFlCgAAAAAsoEwBAAAAgAWUKQAAAACwgDIFAAAAABZQpgAAAADAAsoUAAAAAFhAmQIAAAAACyhTAAAAAGABZQoAAAAALKBMAQAAAIAFlCkAAAAAsIAyBQAAAAAWUKYAAAAAwALKFAAAAABYQJkCAAAAAAsoUwAAAABgAWUKAAAAACygTAEAAACABZQpAAAAALCAMgUAAAAAFlCmAAAAAMACyhQAAAAAWECZAgAAAAALKFMAAAAAYAFlCgAAAAAsoEwBAAAAgAWUKQAAAACwgDIFAAAAABZQpgAAAADAAsoUAAAAAFhAmQIAAAAACyhTAAAAAGABZQoAAAAALKBMAQAAAIAFlCkAAAAAsIAyBQAAAAAWUKYAAAAAwALKFAAAAABYQJkCAAAAAAsoUwAAAABgAWUKAAAAACygTAEAAACABZQpAAAAALCAMgUASGirV6/WoEGDVFVVpfLycg0YMECbNm1yOhYAwAcM0zQbur3BGwEAcIMHHnhAVVVVqqysVGZmpu69916nIwEAvMOIewNlCgCQ6GpqanTOOecoJSVFK1asUCAQcDoSAMA74pYplvnBV2KxmBp5AQFAAtq3b5/KyspUWlqqqqoqp+MAAHyCmSn4hmmaqqqqUiAQUHJysgwj7osMABLMhAkTNHnyZG3btk27du3S3LlznY4EAPCOuBeNQTtTAG6wY8cO9ejRQ8nJyUpKYnIWSHQvv/yykpOTddNNNykajeq8887T0qVLdfHFFzsdDQDgccxMwTdqZ6Zyc3N17rnnSpKSk5MVDPKaAgAAAOLiPVPAdyUlJckwDG3btk3hcJj3UQEAAKDJKFPwLcMwVFhYqEgkopqaGgoVAAAAmoQyBd8zDEOxWExff/21YrGY03EAAACQIChT8D3DMJSUlKTCwkJVV1crEokwSwUAAIBGUabgK4cOHWqwKBmGoXA4zPuoAAAA0CjKFHylrKxMFRUVKisrq/d2wzBkGIZ27Nihmpoalv0BAAAgLsoUfCUjI0MpKSnauHGjioqK6r1P7cYUsVisbtkfAAAAcCzKFHwnEAhoxIgRKikpUWVlpaLRaL33a67t03NyctSnTx+dfvrp+s1vfnNKHtPvpk2bpi5dumjgwIFOR/GUgoICjRkzRv3799eAAQP0zDPPOB0p4VVVVWnEiBEaPHiwBgwYoIceesjpSJ4SjUZ19tln64orrnA6imf06tVLZ511loYMGaLhw4c7HQdwHcoUfCkYDGrQoEEKBALKy8trcNnfqdw+PRqN6s4779SiRYu0efNmvfrqq9q8efNJPSak2267TTk5OU7H8JxgMKgnn3xSmzdvVm5urv7whz/wfD1JLVu21NKlS7VhwwatX79eOTk5ys3NdTqWZzzzzDPq16+f0zE8Z9myZVq/fr3y8/OdjgK4DmUKvmUYhlq0aKGBAwdq48aNCofDDd73VGyfnpeXp9NPP129e/dWixYtNHnyZP3rX/+y/Hg44oILLlCHDh2cjuE56enpGjp0qCQpLS1N/fr1i7s8FifGMAy1bt1akuo2uzEMw+FU3lBYWKh3331X06dPdzoKAB+hTMH30tLSNGLECEUiEX322Wf1vkfqVG2fXlRUpKysrLqvMzMzuThFQti+fbvWrVunkSNHOh0l4UWjUQ0ZMkRdunTRuHHjGNNTZObMmfrtb3+rpCQubU4lwzA0fvx4DRs2TH/605+cjgO4DmccQEeWM7Vq1Urt27dXXl5eg7NPp2L79Msuu8xqVDRgx44dTkfwpLFjx+raa6/V008/rTZt2jgdJ+EFAgGtX79eI0eOVF5enjZt2uR0pIS3cOFCdenSRcOGDdM999zjdBxPWb58udauXas2bdroD3/4gz7++GOnIwGuQpkCviMzM1NnnXWWKioqVFhYWO99Tmb79IyMDBUUFKikpETSkWUpGRkZpyQ7xM6LzSAcDmv16tW6+eabdc011zgdx1MOHTqkMWPG8H6/U+DTTz/V22+/rV69emndunVaunSppkyZ4nQsT6j9HVVaWqqJEycqLy/P4USAu1CmgGOkpaUpFApp//792rhxY72zT1a3Tz/nnHP01Vdfqbq6WjU1NVqwYIEmTJhwqn8E4JQwTVO33367UlJSdNdddzkdxxP27t2rgwcPSpJisZiWLFmivn37Opwq8T3++OMqLCzU9u3b1bt3b1188cWaP3++07ESXnl5uUpLSyUdWZ66ePFidk0FjkGZAuphGIYGDRqkDh06HPXL5FhN3T49GAxq7ty5+uqrr9SvXz9df/31GjBgQHP8CL5y4403atSoUaqqqlJmZqaef/55pyN5wqeffqp58+aptLRUQ4YM0ZAhQ/Tee+85HSuh7dq1S2PGjNGgQYO0ZcsWjRs3jm284Vq7d+9Wdna2Bg8erC+++EI/+MEPWKYOHMNo5OLv1HywDuACpmmqqqpKubm5Ou+88yRJK1asaPS/P/nkEwUCAfXo0UM7d+6Me/9rr71WBw4csPvH8h3DME7ZZ36hYYy1PRhn+zDW9mjfvr3279/vdAzgVIq77WrQzhRAIqr9kN8tW7aosrJSkUhEweDxh86BAwdUUVGhe++9t8E36J5xxhlasmSJevXq1YypG7dv3z5Fo1F16dLF0RxNtW7dOp199tnH/b1bxrU+W7du1ZlnnplwW2CvW7dO1113nT744ANXjuuxSkpKZJqmOnfu7HSURp1++ul14xrvOe1Wpmlq69atrlye+L3vfU8ffvhh3Odroo313r17lZSUpI4dOzqao3fv3lq2bJl69uzZ4P127typkpKShDvXASeDMgWcgGAwqLPOOksfffSR8vLydNZZZ8W971VXXaXHHnss7u0FBQX6/PPPtXz58uaIesLmzJmjffv2afbs2Y7maKpQKFTv2BUWFrpiXOvTrl07LVu2TC1atHA6SpOEQiHNmjVLmzdvduW4Hus3v/mNampq9OCDDzodpVFFRUV14xrvOe1WpmmqdevWrsxcVFSkiRMnxs2WaGP9yCOPKCUlRb/4xS8czVFcXKyrr7660bHLzs62KRHgHrxnCmiC5ORkDRo0SJs2bVJNTQ3LRQAAAHyMMgU0UevWrTVixAhFo1Ft3LjR0nbc7du3b4ZkYFybB+PaPBjX5sG4Ng/GFagfZQqwIBAIqFWrVurcubPy8vIUjUab9P3sNtc8GNfm8Ze//MXpCJ7EuDaPP//5z05H8CTGFagfZQo4Cd27d9egQYNUVVUlSSz7AwAA8BE2oABOUuvWrZWamipJ+uyzzxSNRlVeXh73/jU1NZLU4H3sUFNTo3A47HgOK+rLXF1dHfc2NygvL1c4HHY6RpO5fVy/K5Ge09XV1TJNsy5rImSuVfuikRszV1dXKxaLNZjNjbnjCYfDSkpKcjxzdXV1o7/bAL+iTAGnQO02sJ07d1ZJSYm+/fbbuPetqKg46k+n7N+/XwcPHmwwq1vVl9kt41of0zS1bds2JScnOx2lyUpKSiS5c1yPdeDAAYXD4YR4Th/7fE2EzLVqy5QbM1dUVMg0TVVWVsa9jxtzx3PgwAG1aNHC8cyVlZWKxWLHjevMmTN18ODBuq937Nih4cOHS1Ldn7U6deqknJyc5g8L2IwyBZxC6enp2rFjR4NbpxcUFMgwDGVmZtqY7HjdunVTcnJyg1ndqr7MhYWFkuT4uNbHMAwNHDgw4bZGl6SMjAyZpqmsrCynozSqa9euqqmpSYjndFFRkWKxWN24JkLmWrVlyo2Zi4qKFI1G1aNHj7j3cWPueLp06aKUlBTHMxcXFysSiRw3rh9++OFRX2dnZys/P1+GYSg/P7/Bx+zVq5fS0tIUCAQUDAYbvT/gVpQpAAAA2G7ZsmXq1KmT0zGAk8IGFAAAAABgAWUKAAAAtjIMQ+PHj9ewYcP0pz/9yek4gGUs8wMAAICtli9froyMDO3Zs0fjxo1T3759dcEFFzgdC2gyZqYAAABgq4yMDElHNtmYOHGi8vLyHE4EWEOZAgAAgG3Ky8tVWlpa99+LFy/WwIEDHU4FWMMyPwAAANhm9+7dmjhxoiQpEonopptu0mWXXeZwKsAayhQAAABs07t3b23YsMHpGMApwTI/AAAAALCAMgUAAAAAFlCmAAAAAMACyhQAAAAAWECZAgAAAAALKFMAAAAAYAFlCgAAAAAsoEwBAAAAgAWUKQAAAACwgDIFAAAAABZQpgAAAADAAsoUAAAAAFhAmQIAAAAACyhTAADA06LRqNMRAHgUZQoAAHjWrFmzdMEFF8g0TaejAPAgyhQAoE5OTo4KCwudjuF5W7du1YUXXqjFixc7HcXT8vPz9cQTT+izzz7T22+/7XQcAB4UdDoAADQmEolo27Zteuedd5yO4nm33HKLqqureRW/mW3evFl79+7VxIkTnY7iacnJybriiiu0fv16TZ8+3ek4ADyImSkArjdixAgNHz5cDz74oNNRPK9bt26SpCuvvNLhJN4WjUbVtWtXzZw50+konjZ48GC9+uqr2rJli1auXOl0HAAeRJkC4Hp33323otGorrnmGqejeN6sWbO0dOlS/e1vf3M6iqdNmjRJy5Yt0+zZs52O4hunn3660xEAeBDL/AC43k033aS9e/fquuuu0xtvvOF0HE+bMGGC0xEAAEgYlCkACWHGjBlORwAAADgKy/wAAAAAwALKFAAAAABYQJkCAAAAAAsoUwAAAABgAWUKAAAAACygTAEAAACABZQpAAAAALCAMgUAAAAAFlCmAAAAAMACyhQAAAAAWECZAgAAAAALKFMAAAAAYAFlCgAAAAAsoEwBAAAAgAWUKQAAAACwgDIFAAAAABZQpgAAAADAAsoUAAAAAFhAmQIAAAAACyhTAAAAAGABZQoAAAAALKBMAQAAAIAFlCkAAAAAsIAyBQAAAAAWUKbgK7FYzOkIAAAA8AjKFHxl+/btKi8v1969e2WaptNxHPPZZ5/pzTff1JIlS5Sbm+t0HE+bPXu2IpGIHnnkEaejeNq6dev0zjvv6L333lN+fr7TcTzt0UcflST9+te/djiJt61evVqLFi3S22+/rfXr1zsdB0AclCn4Su/evdWqVSvt3r1bK1euVE1NjS9nq8rKyrR27Vpt3rxZ+/fvdzqOp+Xk5Mg0Tb333ntOR/G0Q4cOacOGDdq0aZMOHjzodBxPe/fdd4/6E81j//792rRpkzZs2KDS0lKn4wCIgzIF30lKStLAgQM1bNgwmaapFStW6JtvvvHVTNWoUaOUmZmptm3b6vLLL3c6jqfNnj1b0n9fzUfzuPDCC9WlSxd17NhRl1xyidNxPO2xxx476k80j/Hjx6t9+/bq2rWrRo8e7XQcAHEEnQ4AOKVly5Zq2bKlRo4cqeLiYpWXl2vz5s3q2bOn09FsMW/ePJWWlsowDKejeNoFF1yg2bNn64orrnA6iqcZhqF58+YpEonwnG5ml1xyiWbPnq0xY8Y4HcXTDMPQ/Pnz1aJFC6ejAGgAZQq+FwgElJWVpYKCAnXq1Emff/65KioqtH//frVv397peM1mxIgRTkfwjZkzZzodwRfOP/98pyP4gmEYPKdtcuGFFzodAUAjWOYHfEeXLl00YsQItWzZUgUFBVq1apXC4bAv31cFAACAhjEzBdQjEAho8ODBqqys1MqVK7Vy5UplZGQ0+L4qwzCUmppqY0r/CoVCTkfwjT59+jgdwRd4TtuHsW5ehmGw1Ba+QpkCGtCqVSulpKRoxIgRKiwsVHl5ub744ot631dlmqYqKiqUm5urkSNHxn3MgoICGYahzMzM5ozeqB07dqiyslJ9+/Z1NEdThUIhlZeXH/f3hYWFkuT4uNbn/fff1/jx4xPuAiMUCmnr1q0yTVNZWVlOx2nUN998o1gspjPOOMPpKI0qKipSLBZTVlZW3Oe0m+Xk5Oiyyy5zOsZxioqKFI1G1aNHj3pvT7Sx/vLLL5WcnKzTTjvN0RzFxcWKRCIVjGgaAAAgAElEQVRxx7VWdna21q5dm3DnOuBksMwPOAHBYFC9evVS69at1bZtW23YsEGVlZU6dOiQ09Ese/PNN/XXv/7V6Ri+MHXq1IS6gEtUr7/+ul599VWnY3hedXW1Jk+e7HQMX5g3b57eeOMNp2MAaAAzU0ATpaenq1u3blq+fLm++eYbRSIRRSIRp2M1mWEYvtoO3klJSUm87w6eEYvFFAgEnI7hC7FYjFkewOUoU4AFhmEoEAho6NChKi8v16pVqyQdWcIXi8VUWVkZ93vD4bAMw2jwPnaIRqMKh8OO57CivszhcDjubU4zDEMVFRVKTk52OkqThcNhmabpynE9VjgcViQSSYisNTU1R41rImSuVV5erqSkJFdmrv0g9oayuTF3POFwWNFo1PHMNTU1rsgBuBFlCjhJoVBIKSkpko784quqqtKXX34Z9/4VFRUyDEOHDx+2K2K9SkpKtH///gazulV9mWt/yTs9rvUxTVNfffWV2rRp43SUJtu9e7dM01RpaanTURq1b98+hcPhhHhOV1ZWHjWuiZC5Vu2SVTdmrh3XsrKyuPdxY+549u3bJ8MwHM8cb1x/9rOfHbXcfceOHRo+fLgk1f1Zq1OnTsrJyWn+sIDNKFPAKdS7d2/t2bNHgwcPjnsft2xA8emnn6qmpqbBrG5VX2Y3b0DRokUL9evXT506dXI6SpNlZmYmzAYU3bp1S5jn9Hc3oJDqf0671aFDhxQMBl2ZubENKKTEGutOnTqpe/fujmeOtwHF0qVLj/o6Oztb+fn5MgxD+fn5jT5uNBrV8OHDlZGRoYULF57SzIBd2IAC8CnDMHgfj00Ya3hJLBZTUhKXD3YwTdPTY/3MM8+oX79+TscATop3j1AADWIDCvuwAQW8hDJlHy9vQFFYWKh3331X06dPdzoKcFI4GwI+lZSURJmyCWMNL4lGo5QpG3l1rGfOnKnf/va3nv354B88gwGfYumZfZiZgpcwM2Ufr85MLVy4UF26dNGwYcOcjgKcNM6GgE8xW2KfpKQkRaNRp2MAp4TX38fjJqZperJMffrpp3r77bfVq1cvTZ48WUuXLtWUKVOcjgVYwtkQ8ClmpuwTCAQYa3gGM1P28Wpxffzxx1VYWKjt27drwYIFuvjiizV//nynYwGWeO8IBXBCmJmyD8v84CWUKft4dZkf4CV8zhTgU8xM2YcyBS+hTNnHqzNT33XRRRfpoosucjoGYJm3j1AADWJmyh6UKXgJZco+zEwB7sfZEPAplvnZhzIFL2FrdPt4dQMKwEs4GwI+xTI/+1Cm4CWxWEyBQMDpGL7gh2V+QKLjCAV8ipkp+1Bc4SUs87MP5w3A/TgbAj7FBb59+JwpeAllyj7MTAHuxxEK+BQzU/bhc6bgJWyKYB/GGnA/yhTgU4ZhUKZsQnGFlzAzZR82oADcj7Mh4FOUKfuwAQW8hN387MMyP8D9OEIBn2K2xD6UKXiJaZrs5mcTlvkB7keZAnyKDSjsQ5mCl7DMzz7MTAHuxxEK+BQzU/ahTMFLKFP24T1TgPtxNgR8ipkp+1Cm4CWUKftQpgD342wI+BQzU/ahTMFLKFP2YZkf4H4coYBPMTNlH8Mw+NBeeAZlyj5sQAG4H2dDwKeYmbIPH9oLL2FrdPswMwW4H0co4FPMTNmHZX7wEmam7MPMFOB+nA0Bn+JDe+1DmYKXUKbswwYUgPsFnQ4AwBmUKfuwpBLHMk1T8+bN0zvvvON0lCZLtDIViUT0+9//XitXrnQ6SpOxzA9wP8oU4FMs87MPM1M41rPPPqs//OEPCXmhnGhlyjRNPfvss9qzZ4/TUZqMZX6A+yXO2RDAKZWIF/iJlrdWIo71gQMHnI7gaVdffbXGjh2bcM8LKfFmS5KTk7VgwYKELCUs8wPcL3HOhgBOqURc5jdnzhynI1iSlJSUMFujm6aphQsXql+/fk5H8bQePXro2Wef1YIFC5yO0mSJOFsycuRI9e7d2+kYTZaIYw34DWUK8KlEK1MHDhzQ73//e0lSVVWVw2maJlFmpkpLSzVy5EjdeOONysrKcjqOL1xxxRVOR2iyaDSqQCDgdIwm27hxo9MRmoyZKcD9KFOATyXapggvvPBC3dKz5cuXO5ymaRLlc6ZSU1N15plnyjRNjRs3zuk4cKlYLJaQZSoRJdqSSsCPOEIBn0q0makZM2Zo9erVkqSzzz7b4TRNkygzU4FAQPPmzdMLL7ygqVOnOh0HLpVoG1AkMpb5Ae7Hbn6ATyXazFQwGFTfvn0lSR07dnQ4TdMk0s6JhmHo+uuvdzoGXIwyZR9mpgD34wgFfCqRLvATXaLMTAEngjJlH2amAPfjbAj4VKIt80tklCl4CRf49mKsAXejTAE+lWjL/BJZIm2NDjSGmSn7sMwPcD+OUMCnWOZnn0AgQHGFZ1Cm7MMsIOB+nA0Bn2Jmyj4s84OXRKNRypRNmJkC3I8jFPApZqbsQ5mCl5imyedM2YSZKcD9KFOAT7EBhX0oU/ASlvnZxzRNyhTgcpwNAZ9imZ99KFPwEsqUfVjmB7gfRyjgUyzzsw9lCl5CmbIPy/wA9+NsCPgUy/zsYxgGW6PDMyhT9mFmCnA/jlDAp5iZsk8gEGCs4RmUKfswMwW4H2dDwKeYmbIPy/zgJVzg24cNKAD3o0wBPsUGFPahTMFL+Jwp+1CmAPfjbAj4FDNT9qFMwUtisRifM2UTZgEB96NMAT7FzJR9GGt4Ce+Zsg8bUADuxxEK+BQbUNgnKSmJ3fzgGZQp+zAzBbgfZ0PAp1jmZx9284OXMFtiH94zBbgfZ0PAp1h6Zh/eMwUvYWbKXow14G4coYBPsczPPpQpeAllyj4s8wPcj7Mh4FPMTNmHMgUvoUzZhyWVgPtxhAI+xcyUfRhreAmfM2UfZqYA9+NsCPgUG1DYh5kpeAmfM2UfNqAA3I8yBfgUy/zsw9bo8BJmS+zDMj/A/ThCAZ9i6Zl92BodXsJ7puxDcQXcj7Mh4EOmaerjjz9WRUWFVq5c6XQcz2OZH7yEMmWPVatWqaysTB9//DHnD8DFgk4HAGC/cDisGTNmqLKyUr/85S/10UcfOR3Jsz744AM99NBDikajOuOMM/SjH/3I6UietHjxYv3xj39ULBbT8OHD9f3vf9/pSJ704osv6rnnnlNSUpL69eunyy67zOlInnXfffeppKREd911lyZPnqyUlBSnIwGoBy8tAT7UokULTZs2TZJ07733OpzG2wYPHqxIJKJoNKpzzz3X6TielZGRoX379unAgQPKzMx0Oo5njRw5UtFoVNFoVEOGDHE6jqfVnpunTJlCkQJcjJkpwKfuuecebdq0SZdeeqnTUTytc+fOGjt2rLZs2aJBgwY5HcezBgwYoDPOOEORSIRxbkb9+/fXaaedpl69eqlbt25Ox/G0Sy65RGPGjOEFL8DlKFPwlZ07d6qyslLffvutUlNT615h9eM2v506ddJ7773ndAxfeP3119nNzwbvv/8+O1TaYM2aNbxnygaGYWjhwoVOxwDQCMoUfKV79+4qKChQq1atVFZWppqaGuXl5SkWi6miokKbN29WKBRSJBJRZWUlSytwSiQlJXHxaYMuXbo4HcEXkpOTnY4AAK5BmYKvBINBBQIBpaenS5L27NmjUaNGyTRNrVixQt26dVN5ebkikYi2bNmiqqoqlZWVad26dQqFQqqpqdGBAwcUCoUc/kkAAADgNMoUoCPLKQzDUIcOHdShQwcVFBRo6NChkqQVK1aoT58+Ki8vV3FxsYqLi1VeXq6ysjLl5uYqNTVV1dXVMgxDqampDv8k/kCZtU+fPn2cjuALPKftw1g3r9rfp4BfUKaAE5CamqrU1FS1aNFCAwYMkHSkZA0fPlwVFRU6cOCATNPUqlWrtHfvXrVt2zbuY1VVVUmS40sIly5dqv3792vSpEmO5miq0aNH65NPPjnu790yrvWZMWOGfve73yXc8qjRo0dryZIlktw5rsdatGiRIpGIrrzySqejNKq6ulqmaSolJSXuc9qtIpGIZs6cqblz5zod5TjfHdf6JNpYv/nmmwqFQho/fryjOeKN689//nMdOnSo7uudO3eqf//+WrNmjYYPH37UfTt16qScnBxb8gJ2okwBJyEYDKpNmzZ1F8lnnXWWcnNz62a16lNQUCDDMBzfvnn58uWKRqMNZnWr+jIXFhZKkuPjWp8vv/xS/fv3V1pamtNRmqxHjx4yTVNZWVlOR2nU4sWLE+Y5XVRUpFgsVjeuiZC5VmVlpbZu3erKzEVFRYpGo+rRo0fc+7gxdzxvvPGGWrRo4Xjm4uJiRSKR48Z12bJlR32dnZ2t/Px8GYah/Pz8uI9XVVWlCy64QNXV1YpEIpo0aZIeeeSRZskONDfKFOBTSUlJisViTsfwhaSkJHbzg2dEo1E2VLFJLBbz5Fi3bNlSS5cuVevWrRUOh5Wdna3LL7+cz+JDQvLeEQrghCQlJbGNtE0CgQBlCp7h14+TcIJXy5RhGGrdurUkKRwOKxwO8z4rJCzvHaEATggzU/ZhrOElXr3Ad6NYLObZkhGNRjVkyBB16dJF48aN08iRI52OBFjC2RDwKcMwuMC3Ccv84CXMTNnHNE3PFtdAIKD169ersLBQeXl52rRpk9ORAEu8eYQCaBTL/OwTCAQYa3gGM1P28XKZqtWuXTuNGTOGnf6QsLx9hAKIi6Vn9mFmCl7CBhT28Wpx3bt3rw4ePCjpyO6QS5YsUd++fR1OBVjDbn6AT1Gm7MMGFPCSWCzGMj+beLVM7dq1S7feequi0ahisZiuv/56XXHFFU7HAiyhTAE+RZmyTyAQYKzhGZQp+3h1A4pBgwZp3bp1TscATgnvvdwB4IRxgW8PlvnBS1jmZx8/vGcKSHQcoYBPsQGFfZiZgpcwM2Ufry7zA7yEIxTwKZb52YeZKXgJM1P28eoyP8BLOBsCPkWZsg9lCl7C50zZhzIFuB9lCvApypR9WOYHL2GZn314zxTgfhyhgE/xnin7MDMFL4lGo8yW2IT3TAHuxxEK+JRhGMyW2ISZKXgJy/zsQ5kC3I8jFPApZqbsQ5mCl5imSZmyCWUKcD+OUMCneM+UfVjmBy9hNz/7mKbJkkrA5TgbAj5FmbIPM1PwEjagsA8bUADuxxEK+BRlyj7MTMFLmJmyD8v8APfjCAV8ig0o7BMIBChT8Aw2oLAPZQpwP45QwKfYgMI+zALCS1jmZx/KFOB+HKGATzEzZR+W+cFL+Jwp+8RiMcYacDnKFOBTzJbYhw0o4CXMTNmHDSgA9+MIBXyKMmUfxhpewgYU9mGZH+B+HKGAT7HMzz5sQAEvYWbKPizzA9yPMgX4FBtQ2IdlfvASypR9KFOA+1GmAJ+iTNmHDSjgJSzzsw/vmQLcjyMU8Cnex2MflvnBS/icKfvwninA/ThCAZ+iTNmHsYaXJOoyv3A47HSEJqNMAe7HEQr4FBf49mGZHxpjmmbCPEcSdZnfokWLJEn33nuvqqurHU5zYljmB7gfRyjgU+zmZ59E3oBiz549KisrczqG540ZM0Z33HGH0zFOSKLOluzcuVOSNHfuXG3YsMHhNCeGDSiO9+CDD+rpp5+u+/r+++/XM88842Ai+F3inQ0BnBKJODO1du1aSUqYV/BrJeLMVCgUkiSNHj1aCxYscDiNd82dO1eStHr1an3yyScOpzkxiToz9cQTT0g6MtYjRoxwOM2JMU2TMnWMadOm6eWXX5Z0pGwuWLBAU6ZMcTgV/CzxzoYATolE3M3v9NNPlyQNGzZMq1evdjjNiUvEmam7775b0pHsffr0cTiNd51//vmSpGuuuUbdu3dPiGMyUTegeOuttyRJffv2dTjJiUvU96c1p169eqljx45at26dFi9erLPPPlsdO3Z0OhZ8jDIF+FQizky1adNGkvTNN9/onnvucTjNiUvEMvXII49IkiZMmKCzzz7b4TTeVTu28+bN00cffZQQsxCJeoE/bNgwpyM0WaIuqWxu06dP11//+le9+OKLmjZtmtNx4HMcoYBPJfJ7pn7605/WvZk8ESTiMr9av/nNb9S6dWunY8BFuMC3D2Ndv4kTJyonJ0erV6/WpZde6nQc+FzQ6QAAnJHIZWr27NlOR2gSPmcKXpKoM1OJiA0o6teiRQuNGTNG7dq147kIx1GmAJ9KxPdMJapEXOYHxJOoG1AkIrZGr18sFlNubq5ef/11p6MALPMD/CoR3zOVqAzDYGYKnpGoG1AkImamjrd582adfvrpuuSSS3TGGWc4HQdgZgrwK2am7MPMFLyE9/HYh7E+Xv/+/fXtt986HQOowxEK+BQzU/ZJ5A0ogGOxzM8+lCnA/ThCAZ+iTNmHDSjgJWxAYR/KFOB+HKGATyXybn6JJhAIsKQSnsHMlH1M0+Q9U4DLcTYEfIr3TNmHZX7wEmam7MNufoD7cYQCPsUyP/uwzA9eQpmyD8v8APfjCAV8ijJlH3bzg5ewzM8+lCnA/ThCAZ+iTNmHZX7wEmam7EOZAtyPIxTwKTagsA9lCl7CzJR9+NBewP04GwI+xQYU9mGZH7wkGo0yM2UTNqAA3I8jFPApZqbsw8wUvISlZ/ZhrAH34wgFfIr3TNmHmSl4Ccv87EOZAtyPIxTwKZb52YfiCi9hAwr7UKYA9+MIBXyKZX724XOm4CXMTNmHczTgfpwNAZ9iZso+LPODlzAzZR82oADcjyMU8CmWntmHDSjgJZQp+1CmAPfjCAV8ijJlH5b5wUtY5mcf3jMFuB9HKOBTlCn7MNbwEj5nyj6UKcD9OEIBn+IC3z4s84OXcIFvH8YacD+OUMCn2M3PPizzg5ewzM8+sVhMhmE4HQNAAzgbAj7Fbn72YazhJWxAYR82oADcjyMU8CmW+dmHmSl4CTNT9mGZH+B+HKGAT82bN0/RaFQvvvii01E8j8+ZgpcwM2WPl156STU1NXr55ZedjgKgAZQpwKeefPJJmaapZ555xukonheJRFRdXa1IJOJ0FOCk8Fy2z7PPPqtoNKqnnnrK6SgAGkCZAnxq1qxZkqSHHnrI4STetnHjRt16663KycnRzJkznY7jWXl5eXrttdf0z3/+UytWrHA6jmfdc889WrhwoaZPn641a9Y4HcfTHn30UUnSww8/7GwQAA2iTAE+ddVVV2natGmaMGGC01E8bcCAAercubMMw9DkyZOdjuNZkUhEX375pb7++mtmTZrRDTfcIMMw1KFDBw0ePNjpOJ52+eWXa9q0abruuuucjgKgAUGnAwBeYpqmTNNUOByOe59oNCrDMBq8j13mzJmjWCyWcO/nqW/sajd4cMO4HutnP/uZnnvuOY0cOdKV+RoSjUYbfU67wTnnnKP09HRFo1GNGjXK9Xmj0ahisVhdTrfnrTV06FD16NFDP/zhD135vIhGo4pGow3mclvmhsyZM8cV43wi4wr4FWUKaCLTNFVTU6NIJKLCwkKVl5eroqJCkrRq1SpVVFRo48aNcb+/qqpKklRSUmJLXi+qb3zdPK7Z2dnKzs5u8HnhVoWFhZKkffv2OZykcY8++qhisVhCjHPt87V2XBMhc6158+ZJcmfm2nHdv39/3Pu4MbfbVVdXyzTN48b17rvv1qFDh+q+3rlzp4YPHy5JdX/W6tSpk3Jycpo/LGAzyhQQRyQSUUVFhcLhsL7++muVl5ervLxcK1euVIsWLRSJRBSNRtWxY0ft3btXknTuuecqNzdXw4YNi/u4BQUFMgxDmZmZdv0onlPf+NZe9DOup1bPnj1lmqaysrKcjtKoho47tykqKlIsFqsb10TK7mZFRUWKRqPq0aNH3Psw1k1XXFysSCRy3LguW7bsqK+zs7OVn58vwzCUn58f9/EKCgp0yy23aPfu3TIMQz/60Y80Y8aMZskONDfKFHyvtiTV1NTo888/V3l5ucrKyrRmzRqlpqYqFospLS1NXbt2VVlZmc477zxJ0ooVK9SzZ09JRz6zyTAMpaamOvmj+EYoFHI6gm/06dPH6Qi+wHPaPox18zIMQ4ZhNHifYDCoJ598UkOHDlVpaamGDRumcePGqX///jalBE4dyhR8Zd++faqqqtKaNWtUXV2tsrIyffXVV0pNTZVhGMrIyFBqaqry8/M1cuRISUdKU9euXSWpwV8QpmmqoqJCubm5dd9bH7fMTBUVFam6ulq9e/d2NEdThUIhlZeXH/f3bp6Z+uSTTzR69GinYzRZKBTS1q1bE2ZmqqCgQLFYrO5FDjf77sxUvOe0m7n1Od3YzFSijfW2bdsUDAYdP/7izUwdKzs7W2vXrm20TKWnpys9PV2SlJaWpn79+qmoqIgyhYTEbn7wlWAwqGAwqP79+2vUqFFq3bq1hgwZojPPPFPJyclq166dWrRo4XRMW7z22mv6y1/+4nQMX5gwYYKqq6udjuF5f/vb3/TSSy85HcPzYrGYLr/8cqdj+MJLL72kBQsWOB2jWW3fvl3r1q1r8EVIwM2YmYKvtG3bVsFgUK1atXI6CgAADTJN0+kIzaqsrEzXXnutnn76abVp08bpOIAlzEwBPub1X9RuwljDK3gu28fLYx0Oh3Xttdfq5ptv1jXXXON0HMAyyhTgU42tacepw1jDa3hO42SYpqnbb79d/fr101133eV0HOCkUKYAAABcyKszU59++qnmzZunpUuXasiQIRoyZIjee+89p2MBlvCeKcDHvPqL2o0Ya3gFz2V7eXEWMDs7m+cRPIOZKcCnvPgL2q0Ya3gNz2l7UDgA96NMAT7GL2oATcV5w14UV8DdKFOAT/EL2l5cgMJLOH8AwBGUKcDHuMC3BxeeAKzgHA24H2UK8Cku8AFYwQW+vThXA+5GmQIAG3ABCi/hAt8enDcA96NMAT7GL2p7cOEJL+G8YS/OH4C7UaYAn+IXNACrOH/Yg+IKuB9lCvAxflHbh7EGYAXFFXA3yhTgU/yCtg9jDS/hhQEA+C/KFADYgAtQeAkvENiD8wbgfpQpAGhmXHgCsIrzB+BulCnAp/gFDcAKZkvsw1gD7keZAnyMX9T2YazhJbwYYx/GGnA3yhTgU/yCtg9jDS/hhQH7MNaA+1GmAABAk/ACgX0Ya8DdKFOAj/Gqp30YawAAvIcyBfgUr3bah7GGl/DCgH0Ya8D9KFOAj/GLGoAVvEBgH8YacDfKFOBT/IK2F8UVXsFz2T6MNeB+lCkAaGYUV3gNz2n7MNaAu1GmAB/jVU/7MNYAAHgPZQrwKV7ttA9jDS/hhQF7cf4A3I0yBfgYF0UArOAC3x6cowH3o0wBPsXFkL24KIJX8FwGgP+iTAFIONXV1U5HaJJELK6HDx92OgJcLBGf04mKsQbcjTIF+FiivcL897//XZL0yiuvOJzEuyKRiJ566in16tXL6SietnnzZv3sZz/ToEGDnI7ieS+88IK+//3v66qrrpIkbdy40eFEJy7RztGAH1GmAJ9KpFc7w+Gwpk6dqttvv12SNH/+fIcTNV2iXBS1bdtWDzzwgGpqapyO4mkbN27UkiVLtG3bNqejNFmiPJdrXX755UpKStLSpUslSd26dXM4UdMk0rka8CPKFOBjiXJRFAgENG7cOL300kuSpNtuu83ZQE2USBdDu3bt0htvvKGxY8c6HcXTJk+erPXr1+uhhx5yOoolifScTk9P1zvvvKM77rhDktSlSxeHE524RDlHA35GmQJ8KpEuhpKSknTLLbfo2muvlSRNnTrV4UTe1aZNG1122WV66623nI7iecFgUHfffbfTMXzBMAzNmTPH6RiWJNK5GvAjyhTgY7zqaR/GGl7Bc9k+jDXgfpQpwKd4tdM+jDW8hue0fRhrwN0oUwBgA15hhlfwXAaA/6JMAT7GRZE9eGUZXsNz2h6cowH3o0wBPsXFEAC4H+dqwN0oU4CP8aqnfRhreAXPZfsw1oD7UaYAn+LVTgBWcf6wD2MNuBtlCvAxXvW0D2MNoKk4bwDuR5kCfIpXO+3DWMNLuMC3F+cPwN0oU4CPcVEEwAou8O3BORpwP8oUANiAiyIAVlBcAXejTAFAM+NiCF7CCwMA8F+UKcDHuCiyD2MNL+EFAntw3gDcjzIF+BQXQ/ZhrOElXODbi/MH4G6UKcDHuCgCYAUX+PbgHA24H2UK8CkuhgDA/ThXA+5GmQJ8jFc97cHFELyE84Z9GGvA/ShTgE9xgW8vLooAWMG5GnA3yhTgY1zg24OLIQBWcI4G3I8yBfgUF/gArOACHwD+izIFADbgAhRewosxAHAEZQrwMS7w7cGFJwArOEcD7keZAnyKC3x7cVEEr+C5bC/O1YC7UaYAH+OiyB5cDMFreE7bg3M04H6UKcCnuBgCYAUX+PbiXA24G2UK8KGioiK99tprWrt2rZYuXep0HF/gAhRewgV+8/voo4+Un5+v119/XQUFBU7HARBH0OkAgJeYpinTNBWNRuPeJxaLyTCMBu/T3EpLS7V48WJJ0ueff64LL7zQsSxW1Dd2sVgs7m1Oq/3/243ZGhOLxRp9TrtFLBZTLBZLyKyJkLmWm4+1E3kOuDF3fbZs2aKtW7dq69atOnz4sKO5E+nYAuxGmQJOoVWrVqmiokJr1qyJe5+amhpJ0n/+8x+7YtXrzDPP1LfffqshQ4Y0mNeN6svrlnE91qFDh1RRUaGlS5dq4MCBTsdpsp07d0qSdu/e7XCSxhUXFysSiSTE87n2+Vo7romQudYXX3yhqqoqffjhh2rXrp3TcY5SU1Mj0zS1Z8+euPdJlLEeOHCgkpOTlZWVpfLyckdzxxvXe+65R4cOHar7uqCgQMOHD5ekuj9rderUSTk5Oc0fFrAZZQo4CTU1Ndq5c6fKysokSUOHDtXatWs1YsSIuN9TUFAgwzCUmZlpV8x6zZ8/X2vXrtXo0aMdzWFFfeNbWFgoSY6P67HOP/987d+/XzNnztTu3bvVqlUrpwmiciUAACAASURBVCM1Sa9evWSaprKyspyO0qDly5dr2bJlisVimjx5si666CKnIzWoqKhIsVisblwbOme4STgc1rhx4xSJRPTAAw8oLy/P6UhHKSoqUjQaVY8ePeLeJ1HGWpKee+45DRgwQIMHD3Y0R+0LFceO67///e+jvs7OzlZ+fr4Mw1B+fn6Djzlt2jQtXLhQXbp00aZNm051ZMA2vGcKsKC8vFyVlZXKz89Xy5YtFQqFJEktWrRwONmJGzBggKZOnep0DM+78847JUmXX355whWpRJKSkqIdO3aooKBALVu2dDqOZyUnJ+uqq66SJP3kJz9xOI333XTTTY4XqeZy2223MVMFT6BMASfINE1FIhGtXbtWmzdvVnJyskaNGqWsrCzejI24rr/+emVlZelXv/qV01E8bfjw4crKylJ6erpGjRrldBxPu//++5WRkaGbb77Z6ShIYBdccIE6dOjgdAzgpLHMDzgBu3bt0o4dOxQOh/W9731Pbdu21YoVKyhRaFQwGNQXX3zhdAxf+Ne//qVIJOJ0DM/r06ePvvzyS6djAIArUKaAOCKRiAoKClRWVqbDhw9ryJAhWrt2rdq2bet0NAD1OPPMM52OAADwGcoUcIzKykpVVVVp1apVysjIUCgUUp8+fRr9PsMwlJqaakNC1L5HDc3vRJ77OHk8p+3DWDcvwzBYtQFfoUwB/79oNKr169erqqpKgUBAo0aNUlJSkoqLi0/o+03TVEVFhXJzczVy5Mi493PLbn4vvfSSQqGQJk2a5GiOpgqFQiovLz/u7926m191dbWmTp2q1157zekoTRYKhbR169aE2M1Pkv7+978rHA5rypQpTkdp1Hd384v3nHazG2+8Uc8//7zrXkBqbDe/RBvrt956S/v27dPtt9/uaI54u/kdKzs7W2vXrqVMwVfYgAK+t2fPHuXl5am6ulo9e/bUyJEjlZycrKQkbx8eW7ZsUVFRkdMxPO/QoUNatWqV0zF8Ydu2bfr222+djuELq1evPurzhdA8ioqKtGXLFqdjNIsbb7xRo0aN0tatW5WZmannn3/e6UiAJcxMwZei0aiKiopUVlamkpISDRgwQBs2bFD79u2djmabmpoapaSkOB3j/2vv3qOjqu/9/79mksxMLhBCuF8CuZAAllsgKFKkWiyISKm3Sm27TqueWrX21EtL+y12HXQV6zq1arVVW4+iFqyKFloppwUFKyqIchEJSSAh5EK4JQTI3DIz+/cHv6SiBGGEvffMfj7WcrlgZoV3NiGZ17zf+/1JeoFAgJXoSDrp6eny+/1Wl5H0fD6fQqGQ1WWcE0uWLLG6BOCsSO633oFPCIVCCoVCeuedd9Te3q7MzEyNHDnSkTP0wWCQ83hM4Pf7bTcKBXxehClzeDweBYNBq8sAcAp0puAo+/btk8vl0oUXXii32619+/ZZXZJlQqEQYcoEwWCQDiCSTnp6Oi/yTeDz+RQOh60uA8Ap0JmCo+Tl5cnj8ST9/VCngzBlDjpTSEYZGRl0pkzg9XoJrYDN8YoScCg6Jubgnikko/T0dAUCAavLSHqEKcD+CFOAQ4VCIXk8HqvLSHqEKSQjwpQ5GPMD7I8wBThUKBSiM2UCv99PmELSYQGFOehMAfZHmAIcijE/cwSDQe6ZQtLJyMigM2UCwhRgf4QpwKEY8zMHnSkkI5/PR5gyAWN+gP0RpgCHYszPHIQpJCM6U+agMwXYH2EKcChWo5sjGAwSppB0WEBhDp/Pp1AoZHUZAE6BMAU4VDAYJEyZgM4UkhELKMzh8XgIU4DNEaYAh2LMzxwsoEAyYszPHD6fjzE/wOYIU4BDMeZnDjpTSEYsoDAHY36A/RGmAAcyDIMwZRLCFJIRnSlzpKamyjAMRSIRq0sB0AXCFOBA7e3tSk1NldvNt4BzjQUUSEbcM2UeulOAvfFKCnAgDuw1D50pJCO2+ZmH9eiAvRGmAAdik595WECBZMSYn3m8Xi+dKcDGCFOAA4XDYcKUSehMIRn5fD7G/EzCmB9gb4QpwIEY8zNPIBAgTCHpZGRkMHpmEsb8AHsjTAEOFAwG5fF4rC7DEQhTSEYsoDAPYQqwN8IU4EDhcJjOlEkY80MySk9P5wW+SXw+n8LhsNVlAOgCYQpwIBZQmIcFFEhGGRkZdKZMQmcKsDfCFOBAHNhrDsMw6EwhKXm9XoXDYUWjUatLSXps8wPsjTAFOFAoFGLMzwThcFgpKSlKTU21uhTgrHK5XIz6mYQwBdgbYQpwIMb8zMHyCSQzllCYw+fzEVoBGyNMAQ7EmJ85CFNIZunp6bY7uHfOnDn68pe/rDvuuEP19fVWl3NW0JkC7I0wBTgQnSlzBAIBlk8gaWVkZNguTF100UVqbGzUE088oYaGBqvLOSsIU4C9EaYAB2I1ujn8fj/X2YGCwaBef/11tba2Wl3KOeXz+Ww35nfbbbcpNzdXEyZM0E033aRNmzZZXdLn5vP5CFOAjRGmAAfi0F5z0JlynlmzZqlPnz6aPXu2/vCHP6ipqUmGYVhd1jlhx86Ux+PRqlWrtHbtWs2fP19z5szR448/ntB/Bx6Ph3umABtjxRTgQMFgkI6JCQhTzjN37ly9++67CgaDWr58uR5++GFJ0vDhwzV8+HCNGDFCw4cPV05OjnJzcy2u9vOx4z1Tkjq/t11zzTUqLS3Vt7/9bQUCAf3oRz+yuLL40JkC7I0wBThQOBzmnikTJNKYn2EYampq0s6dO1VWVpYwddvN9ddfr6985St65plndOedd8rlcmn//v3asWOHysvLtWPHDi1btkzbtm3TtGnTtGjRIqtLjlsibPMrLCzUmjVrEjqM+Hw+tbW1WV0GgC4QpgAHCgaD6t69u9VlnDbDMLRixQpdeumlVpdyRuzUmfL7/brjjjs0btw4zZo1q/OFfcf/d+zYodTUVI0cOVJPPvmkBg8ebHXJCat37966++67O3/dt29f9e3bV1OnTu38vYaGBkUiESvKO2vsOOZ3MmlpaUpLS7O6jLh5PB41NzdbXQaALhCmAAdKpEN7m5ubdfPNN2vv3r26+OKLrS7njAQCAVtc58rKSn3xi19UW1ubFi9erIULF2rEiBEaMWKExo4dq+uuu07Dhw9X7969rS7VUdzuxL5t2efzJUSYSnQ+n0/hcNjqMgB0gTAFOJDdw1RbW5u2bdsmwzD0H//xH/rqV7+q559/PuGWZtilM5WVlaUJEyZo69at6tWrlzZv3mx1SUgCidKZSnQc2gvYG2EKcCC7b/NbuHChHnroIWVnZ+vJJ5/U5ZdfbnVJcbHLob0DBgzQihUrZBgG73DjrLHrAopkwzY/wN4Se8YAQFzs3JlqbGzUww8/LMMw1LNnz4QNUpJ9wlQHl8vF4hGcNYQpczDmB9gbnSnAgUKhkG1fVPv9fk2ZMkXXXHONLrzwQqvL+Vz8fr8txvyAcyE9PT3pDya2A6/XS2cKsDHCFOBAdj5nqqioSCtWrLC6jLMiEAgoOzvb6jKAcyI9PV1NTU1Wl5H0CFOAvTHmBziQnTtTycQuCyiAcyEjI8P250wlA8b8AHsjTAEORJgyh93umQLOJu6ZMgedKcDeCFOAA9l5AUUyIUwhmRGmzOH1ehUKhawuA0AXCFOAA9l9NXqyYAEFkhlhyhyEKcDeCFOAA9GZMgedKSQzDu01B4f2AvZGmAIciDBlDsIUkll6ejoLKExAZwqwN8IU4EChUIgxPxP4/X7CFJJWeno6HRMTEKYAeyNMAQ5k53OmkkkwGCRMIWnRmTIHY36AvRGmAAdizM8cLKBAMmMBhTm8Xq/C4bAMw7C6FAAnQZgCHCYWiykSiSgtLc3qUpIenSkkMxZQmMPtdis1NZWDewGbIkwBDtNxYK/L5bK6lKTHPVNIZoz5mcfn83HfFGBThCnAYYLBoLxer9VlJL1YLMY4JZJaWlqaXC6X2tvbrS4l6Xm9Xu6bAmyKMAU4DGHKHB1LPtxuvs0ieWVkZNCdMgEb/QD74qc84DDhcJhuiQkY8YMT+Hw+7psyAWN+gH0RpgCHoTNlDg7shROwhMIcHo+HMT/ApghTgMMQpsxBmIITsB7dHHSmAPsiTAEOw5ifORjzgxMQpsxBmALsizAFOAydKXMEg0EO7EXSYwGFORjzA+yLMAU4TMc5Uzi3/H4/HUAkPRZQmMPn83FoL2BTqVYXACSTTZs26dixY3r33Xe7fE44HJbL5VJ9fb2Jlf3b1q1b5ff7T1mj3Z2s9o6zbqy6rp+0efNmhUKhhL7OklRTUyPDMNTQ0GB1KZ+prq5OkUgkIa55e3v7Cdc1EWo+mWAwqC1btqh79+5WlyLp39e1sbGxy+ck4rVua2vT1q1b1aNHD0v+/K6u67x589Ta2tr567q6Ok2YMEGSOv/foVevXlq5cuW5LxYwGWEKOAui0agkaeDAgQqFQjr//PO7fG5dXZ1cLpcGDRpkVnknaGhoUP/+/XXBBRdY8uefDServSNEWXVdP6m+vl4DBgxI6OssSfn5+TIMQ4MHD7a6lM+0Zs0ahcPhhLjmDQ0NisVindc1EWo+mQEDBmjw4MG2qb+hoUHRaFR5eXldPscutZ6J/v37Ky8vz7LaGxsbFYlEPnVd16xZc8Kvv/jFL2rjxo1yuVzauHHjKT/mypUr9cMf/lDRaFQ33nij5s2bd7bLBkzBmB/wORiGoerq6s5Z9j59+lhc0am1trbq7bff1sGDB7V3716ry0lqLKCAE7Aa/dxramrSwYMH9e677+rw4cNWl3NWRKNR3Xrrrfr73/+u7du3a8mSJdq+fbvVZQFxIUwBcTIMQx988IHC4bAyMzOtLue0fPjhh3r88cf11ltv6ZlnnrG6nKTGAgpztbe369ChQ2pububeEhOlp6ezgOIcW7x4sdasWaM//OEP2rRpk9XlnBUbNmxQUVGRCgoK5PF4dN1112nZsmVWlwXEhTAFxKGlpUVtbW0aPHiwhg8fbnU5p23y5Mnq2bOnUlJSdPPNN1tdTtJauXKl5s+frxdeeEFLly61uhxHWLJkiX7/+9/rj3/8o5599lmry3GEZcuW6fnnn9eCBQv0t7/9zepyktYNN9yglJQU9ejRQ1OnTrW6nLOioaHhhNHhQYMGJcR9mcDJcM8UcIaqq6t14MABZWRk2H6s75NcLpd+9rOfac+ePcrJybG6nKQ1YsSIznfrhw0bZnE1zjBnzhzdfvvtkqSrr77a4mqcoaioSG1tbZKOf83j3MjOztbtt9+u3r17y+3mPXDAbghTwGkKh8Py+/1qb29XWVlZlxuhDMOQYRif+fFO5znnQkdHyqo//2w4We0dv2eHzysvL0+FhYXy+XwaNWqULWr6vOz+OXTr1k3XXXedwuGwsrOzbV9vx/cJO33dnqmRI0dqzJgxam1tVUFBga0+h1PVYqc6T9eCBQskWVf7J79eTyYUCp12fQMHDlRdXV3nr+vr6zVw4MDPXSdgBcIUcBqam5tVXl4uj8ejkpKSLp/X8cOmq8MVI5GI9u7dq0GDBnEA4+dwsmvX0tKinj172ua6rly5UmlpabapJ16BQCBhPo/f/OY3kk7+9WE3sVhMwWCws9ZEqPlkli1bpnA4bKv6W1pa1Lt37y4ft1OticLtdqupqUk5OTlKS0s76XNWrFihuro6HThw4DM/XllZmaqqqlRTU6OBAwfqhRde0OLFi8922YApCFPAKRiGoVAopKqqKpWWln7mzb/RaFSpqanasmXLSR8LBoPyer3as2eP9uzZc67KTnqfvL7t7e2KRCLy+/22OWcqWdTX18vr9aqpqcnqUpJKR5hqbm6W9Omv6URTW1trdQmdAoGANm7c2OWL/kS/1laJRCJ677335PP5lJKS0vn7Hz9rKjs7W0OGDJH06XOmpH+fNZWamqpHH31U06dPVzQa1Xe/+12dd9555nwiwFlGmAK6EA6HtXXrVhmGobKysi5n1Tu6UT169LDNwZXJzOVy6ctf/rLVZTiCy+XSrFmzrC4j6fE1bR6utTlSUlJ0ww036Oabb5bL5Trpc2bOnKmZM2eaXBlw9hGmgJPoeAeuuLhYVVVVpwxSkUhEKSkpXW5tC4VCikajSk9P7/KHihUqKys1dOhQeTweq0uJm2EYCgQC8nq9J7xTaifhcFi7d+9WcXGx1aWcsWg0qnA4nDDnZR08eFCGYZxyxMtOAoGAPB6Pbb92T6WyslJDhgyR1+u1upST6pgEyMjIsNX33TMVDodVXV1tq62xHd93U1JSPvX339GlMgxDCxYs0Pz585WXl/epn6EdHSogGRCmgI/pOIQ3FAppypQp8vl8qqqq6vL57e3tkqTS0tJPPRaJRFReXq6cnBwVFhba7gf6bbfdpqeeekpf+MIXrC4lbjt37pTH41FeXp7VpXRp5cqVeuWVVxJyRXpdXZ3cbnfC3Bj+wAMPKBwO6+c//7nVpZyWxsZGRSIRW3/9duXXv/61rr/+el1xxRVWl9Kluro6BQKBhHwjo0NlZaV+8IMf6IMPPrC6lBN0/Kxsa2vTyJEjlZp6/OXk6tWrO58Ti8X01FNP6c9//rMWL16swsJCq8oFzinCFPD/MwxD77//vrp3767MzEz5fL4unxeLxeRyufTee++d9DmxWKzzXedgMKiNGzeey9LPWDQa1a5du9Tc3Nzl52B3kUhE4XBYGRkZ2rdvn9XldOnNN99Uenp6Ql7njq5fY2Oj1aWclvr6+s6uciKIxWIKhUK2/vrtitfr1b/+9S/169fP6lJOKRAI6NChQ50v9hNNOBxWbW2t3nnnHVt+Du3t7XrrrbeUnp5+Qvfp4/dRxWIxnX/++crJyVF2dvanPgZdKiQ6+/3LBCxw6NAhtbW1qbi4WL1799bBgwe7fG40GlUsFtP48eNP+nhzc7N27typ0tJS295DVVNToz59+uiiiy6yupS4hEIhbdmyRZMmTbL9mOIrr7yisrIylZWVWV3KGTEMQxs2bNDEiRNt11Xtyuuvv65wOJww17rjGk+YMCFhrnGHiRMnqqamxvbXur29XZs2bdKoUaO6fIPM7gYOHKhevXrZ9sy6o0ePqry8XAUFBcrNzZV0YodKkg4cOKAbbrhBU6ZM0S9+8YuEHG0FukKYguPt3LlTzc3NysjIOOW9FrFYTG63W+vXr+/yOeFwWJFIRD6fTxUVFeei3LNiw4YN6tOnT8K8g/9Jfr9fHo8nIbZybdq0SZdeemnCXeuOronduqqn0tDQkFCdKen4GwPr169PuBeXkUhEW7ZsSYhrHYvFtH79etvdt3q6+vTpoxUrVujCCy+0upRT2rZtm1JTUz/1BtfHu1TPP/+8nnjiCfXv3/+knTa6VEhEhCk4VigUUltbmwzD0IQJE7o8hFc6/sM4Eomc9N6ojsc7wlNJSYntT6l/7733NH78eNu/q3wy1dXV6tWrl/Lz860u5bS0trZq2rRpGjNmjNWlnJGOFfODBg2yuJLTl2idKel4AIzFYho8eLDVpZyRzMxMPffccwlzrXfv3q1oNJqQ9+2MHz9eLpfL9tc6FoupqqpK0WhUJSUlnW8QfLxLZRiG/va3v+n+++/X7373O02cONGqcoGzhjAFRzp06JB27Nghr9d7ytEJwzDkdrv1zjvvdPmcjvNiOt6Re//9989FyWfV22+/rcGDByfEu8of9/H7pE41imkXHTdpJ+K9aR33/O3du9fqUk5bInamOjqAiXaOVyAQ0O7du7V+/Xrbv3nUwe/3a//+/ba89+hUPB6P3n77bU2aNMnqUk5LOBzWW2+9JZ/P96mvjY4uVSwW0zXXXCO32905GvhJdKmQKBLrOwrwOXUcwltdXa3x48efckNSx9rzrrpR0vGuQ0VFhUaPHq2cnJxzUfI5ceTIEX35y1+2/TudHxcOh7V582ZdcMEFtl3H/EkHDx6U1+vVxRdfbHUpZyQR75eSErMzlcj3TWVnZysvL0/9+/e3upTT0nGv5ZgxY2x/r+XH+f1+bdy4MaG+rg8fPqzKykoVFxefsHTi412qY8eO6fbbb1f37t312GOPJew9bUBivJ0EnCWVlZWSjp/M/lnfuNvb22UYRpeP7927V1VVVRo1alRCBSlJqqqqsu3NzCdjGIbKy8tVWFiYMEFKOj5alCjjiB8XCAQS/nyeROFyuZSZmSm/3291KWds6NChqq2ttbqM0+b1elVUVKTt27ef8nu73RQVFZ3yiA476tGjh0aPHq2qqqout4FmZWXpj3/8o0pKSnTxxRerpqbG5CqBs4POFBxl2LBhamlp6fJF4umsPe/obhmGIZ/Pp23btp3Lks+6QCCgAwcOaO/evdq/f7/V5ZyWcDjcOTJXXV1tdTmn7fXXX1e3bt0SauxMOn69JSVc3Yk45if9u+uaSN0SSerWrZtef/31hFueEQqFtG7duoR5Y8YwDB05ckRvvPGGsrKyrC7nU44cOaJ7771X+/btU9++fXXPPfeoW7duko7XvnPnTu3atUtr167Vn/70J0nS9ddfrzfeeEOtra0KBoNqampSUVGRunXrpqKiIklSS0uLGhsbFQwGdcEFF5wwbr9w4UI99dRTSklJ0SOPPKLp06dLOn6u3w9/+ENFo1HdeOONmjdvnslXA05EmIKjfNZs/2etPW9vb9dHH32kPn36aMiQIQn5zv2WLVtUVFSkCy64wOpSTktra6t27typcePGJcy9GR3Wrl2r0tLShBrPkaSPPvpIeXl5nS+IEkUijvlJx8eddu/enXAHaJeWlibEYoRPisVi2rx5swoKCtSjRw+ryzktJSUl6t69e5c/m6w0f/58zZ49W3fccYcefPBBrVmzRgsWLOh83DAMffjhh3ruuef05ptvyuPxaOrUqVq7dq1ycnL0pS99Sf/7v/+rvLw8nX/++Ro+fLgWLVqkyspKud1ufe9739P//M//dH687du364UXXtBHH32kxsZGTZs2rXPq5NZbb9U///lPDRo0SGVlZZo9e7ZGjhxp+jWBsxCmAJ3e2vNoNKpgMCiv16v29nYdOHDAxArPnjfeeEO5ubkJ8e69YRjy+/1KT09PiMUen7Rx40YVFxcnxLX+uLa2toQcO0vUzpR0/JonYt3vv/9+QtbdEagSZZy1Z8+e+r//+z/FYjGrS/mUpUuX6sEHH9R7772nvLw8LVy4UC+99NIJXarXX39d48aNU3l5uXw+n7KzszV27Fh5vV6Fw2H99Kc/VWtrq9xut5YsWaIXXnhBOTk5GjJkiCoqKnT99ddr7969CofDSklJ0Z133imv16uXX35ZBw4c0IgRIxSJRFRTU6MePXrI4/Foz549uuSSS9SvXz+lpqYm1DEPSCyEKTjeZ609l44fOLh7926VlZUpMzPTxOrOvtWrV+v888+3/bvJhmFo27ZtGjp0qPr06WN1OXFpa2vT1KlTbX+tPy4QCKiqqkqjR4+2upQzlqidKUn68MMPVVBQkFDfXwKBgN5+++2EvN6SOsedR40aZftANXHiREWjUVte66NHj2rGjBmSpGXLlkmSysvLT+hSvfPOOyotLdXEiRP1zjvvqKamRrfddpsmTZqka665Ri+//LJycnI0YcIEFRUV6fvf/75+8IMf6Nvf/rZeeeUV5eXlaeTIkZo3b54mTZrUuTzq7rvv1o4dO3TZZZdp48aNevbZZ9WzZ09JUmpqqi6//HI99dRT1lwYOAZhCo5lGEbnxr5TPaempkZHjx7V2LFjlZaWZmKF50ZVVZUuueQSq8v4TA0NDfJ4PAkbpKTjCyiGDh1qdRlnpKWlJWFGn5JJjx49dPjw4YQKU0OHDtXu3butLiNuvXv3VktLi+rr621/ztewYcO0YsUKy/782bNna9++fZ/6/XvuueeEX69YsaLzXrSZM2dq6tSpWr58uVwul+bMmdN5rMWAAQMUCASUlZUlr9ersrIypaen68iRI8rPz9eVV16pd999V3fccYfcbrfef/993X///ZKk/Px8LV26VMOGDdPQoUPVt29fSdJrr72m/fv3a+zYsZJ0Qnede6lwLhGm4EgdISolJeWUiyaCwaDcbre8Xq82b95scpXnxtatW3XRRRfZejSnY6QyIyPD1nWeSiQSUWNjo/bt26fm5maryzltHedLneyFk90l8phfNBpVOBzucvOZHUWjUR04cEDr1q1LuOUZHQzDUGNjoxoaGmy9SKO9vV1btmyx7Gt7/vz5XT7WrVs3rVy5Urm5uWpsbFTPnj313nvv6YknnlAkEtEf/vAH3XfffXr11Vc1c+ZMbdy4UbFYTGlpaSovL9eRI0f07W9/W6NGjdK8efO0ceNGTZgwQbt371ZaWpqi0aiOHTumiy66SD169FBFRYXC4bCys7NVXl6ut956SzfccIN27dqlHj16dP6s7tmzp1auXKnS0lLV1dVp/fr13EuFc4IwBUdqb2+XpC5H+/x+vz766CMVFxerX79+ZpZ2TnW8cLjiiis6RyHsJhKJaNOmTQk/UllTU6N+/frpwgsvtLqUM9Jx5lGiLfuQEnvML1HPmxo0aJD69OmTUEctfFLH9/tx48bZ9kDfkpIS3XnnnRo/frzl/zY/2aXy+/269dZb9cgjjygajeqqq65SWVmZvve97yk9PV1lZWX61a9+pQsvvFBFRUXKzc3Vvn37dMstt+jVV1+Vx+NRSkqKxo8fr5SUFH3zm9/UrbfeqksuuUTnnXeejh07pnfffVfHjh3T6tWr9aUvfUl+v1/r1q3T5s2bNWXKFDU1NWnMmDHasmWLampqNHDgQPXr108vvfSSdu/ereuuu0719fUqKCjQddddp2XLlhGmcNa4PuOshcQ5iAH4DIZhqK2tTRs3blQ0Gj3lc0OhkNLS0iz/oXW2HTp0SP/5n/+ppUuXWl1Kl6LRqAzDsO2LmtP1wQcf6E9/+pN+/etfW13KGWlvb0/Ycdbn4NCguQAAFydJREFUn39e7e3t+s53vmN1KXFJxGv/k5/8RFdddZUmTpxodSmfSyQSkcvlsnV36utf/7p++9vf2m70ubW1Vffdd5/279+vgwcP6rHHHtPQoUM1a9YspaSkaNmyZWptbdXVV18tt9stl8ul4uJiPfLII3rxxRe1cuVK+f1+ud1uHTp0SJmZmfJ6vQoGg2pra+scyZeOb+Q1DENut7vz7yocDsvtdis1NVXRaFQpKSkyDEPdu3dX//79FQ6H5XK5dOONN+quu+5S7969FYvFNHjwYBZT4Ex0+S5XYr9aAc5Ax+GYU6dOtboUy6xZs0Zf+MIXEuKeqUS3a9cujRs3jmttonXr1ikcDnPNTVRaWqru3btzzU0watQo9erVy1bX+tJLL1VTU5Ok4/f9BYNB/exnP9MTTzwhSfr+97+vSy65RD/+8Y/l9Xrl9/v1i1/8Qg8//LDGjRundevW6ciRI9q5c6ckqW/fvrr22mv19a9/XdOnT9fixYt17bXXqlevXsrKylJlZaUyMjI0Z84cvfzyy/rxj3+sRx55RBkZGbr55psVCAT0m9/8Rr/85S/15JNP6sMPP9Tdd9+thx566ISjB+bMmcNiCpw1hCk4SiKNz5wLVVVVKikpcfx1MENNTY0KCwu51iZyuVyd/8EchYWFqqmp4ZqboLi4WJWVlbr00kutLqXTqlWrTvj1oUOHdO211+ruu++W2+3Wd77zHblcLi1ZskQul0vDhw9XOBxWLBbTxIkTdfToUY0ZM0a5ubl69tln5fV69fTTT+vZZ59VSkqKotGompqaFA6H1dDQoG7dusntdndOjTz33HOKRCKdo/srVqzQb37zG02ePFk///nP5fF4lJqaqlgsphkzZigUCqm1tVWvvvqqtm3bpj//+c8JtyQI9pNcM0wATqmiokLFxcVWl+EINTU1KigosLoM4JwqKChQTU2N1WU4QklJSefhtHaVm5ur1atXq6qqSrfccouWL18uSTp48KA8Ho/+/ve/a8eOHfL7/Vq+fLnKy8tVW1urlpYWvfjii4pEItq/f78WLVqk1NRUvfPOO3rzzTcVCAS0dOlSVVRUyOv1asuWLYpGo9q3b5+++93v6tChQ1q2bFnnApd7771XPXv2VCgU0ujRo5WWlqaamho98cQTikaj6tevn/bv36+rr77aysuFJEFnCnCQiooKTZkyxeoyHKG6ulr5+flWlwGcU/n5+aqurra6DEcoKSnRP/7xD6vL6NK0adM6R/6k4/e/1tXV6eGHH+48y7GgoEAbN26Uy+XSsmXL9NOf/lTz589XWVmZ9uzZo0mTJqlnz566+uqr9Z3vfEcvvviinn32WWVmZmr27Nmd91vt2rVLQ4YMkSTdf//98ng8KiwsVHl5uUaPHq3y8nKtXr1aa9asUVNTk0pKSjR9+nTV19dr7ty5ev7559XY2Ki8vDytXbvW0eP/+PzoTAEOUllZqZKSEqvLcITq6mo6U0h6BQUFqq6u1mcss8JZUFxcrIqKCqvL6NKqVau0bdu2zv/Ky8t17Ngx7d27V3369OncjDtw4EB1795dDzzwgIqKirR//37t3LlTffv21aBBgyQdD2KRSESGYcjj8cjj8XSOkh4+fFjt7e2KxWJKT09XKBTq/Lgul0tHjx5V9+7d1aNHD914440aNmyYduzYIY/HI6/Xq3vvvVfS8Z+H0WhUX//61zV27FgtWLDAgquGZECYAhwiHA5rz549KiwstLqUpHfkyBEFAgHbbd0CzracnBxJxw97xrmVn5+vvXv3KhgMWl3KGZswYYKqqqokSU8//bRCoZBmzZql7du3a8mSJdq+fbuysrL05ptvqqWlRb/97W8VjUa1bt06lZaWKhAI6N1339VHH32kuro63Xbbbfp//+//KRKJ6Omnn1Y0GtVzzz2n4cOHa926dTp69KhmzZql5557Tt/85jdVVFSkDz/8ULFYTM8884yk42fqud1uPfzww9q8efOnDiAGThdhCnCI6upqDRo0KGEP10wkNTU1ys/P56Z8JD2Xy9XZncK5lZqaqvz8/M7Nd4nklltuUXV1tYYNG6ZXX31V/fr108iRI7V161YFAgFNnTpVtbW1KioqUllZme655x75fD4VFRVp0qRJCgQCuuiiizRmzBilp6drzpw5Ki0tldvt1oIFC+Tz+dTe3q4rr7xS/fv3VzQaVX19vS6++GLdeeedqq+vl8vlUiwW03333af09HR99atflWEYuvbaa2UYhm6//XYVFRVp9OjR+uCDD6y+ZEgghCnAIRjxMw8jfnASwpR57D7q15Vp06YpOztb//jHP3TnnXeqpaVFs2fP1rhx49Ta2qoZM2Zo4cKF2rp1q5YvX67c3FxdcsklcrlcnedC1tTUaNKkSQqHwzr//POVmZmpcDis+++/X7/+9a8ViUR0zTXXKBaLKS0tTZMnT9bRo0c7g1LHaGBHp2vevHlKS0vT2LFjVVZWpk2bNqmqqkpPPvmkvv/971t8xZBIWEABOERFRQVhyiQdnSnACfLz89noZ5JE2Oh3MqmpqXr00Uc1ffp0tba2qqCgQOedd56++93vyufzKTs7W9/73vf0xz/+UVOmTFEgENBPf/pTSdKGDRuUmZmpyZMnq6WlpXM1+muvvaaCggL96le/UiwWU1ZWltra2rRhwwYNGTJElZWVGjlypAYMGCCPx6NbbrlFeXl5Sk1NVVFRkbKzs7Vx40aNGjVKl112mTZt2iSXy6ULLrhAhw8f1t69e9W/f3+LrxwSAZ0pwCFYi24eOlNwEjpT5ikpKUnIzpQkzZw5U5WVlfrLX/6i3Nzczt/r06ePBg4cKJ/PpxkzZigWi0mSnn/+eUlSQ0ODUlJSVFNTo7vuukvRaFQFBQV66KGH5Ha79corr+gvf/mLJOlrX/uarrrqKhmGof/6r//Szp07NWPGDP3yl7/UkiVLlJ6erubmZmVkZKipqUnf/OY3NXbsWL3xxhs6fPiwDh48qDVr1mjXrl2aOnUqiylwWghTgEPQmTIPZ0zBSThryjyJOub3cWVlZaqqqlJNTY0ikYh27dql2bNnKxqN6tlnn9WsWbO0cOFCvfjii9q+fbtaW1s1efJkuVwuhcNhRSIRvfXWWyouLlZtba0mTpyozMxMHTt2TPPmzdO8efO0a9cuXX755YpGo3riiSeUmZmphoYGhUIhvfLKK3rvvfc0efJkXXnlldq8ebOysrLk8Xg6Q16PHj20ePFiFlPgtDDmBzgE90yZhzOm4CScNWWejs6UYRgJu+Dm4yN/bW1tys3N7Rz5S0lJ0ciRIztH/s4//3z5/X61t7fr/vvv1/r165WWlqYpU6boyJEjMgxD3/jGN7Rq1SqlpKToV7/6ldxut7xer37/+9/rr3/9q44ePaqsrCxNnTpV9fX1Gjt2rCQpFotp6dKleumll9Tc3Kzbb79dX/nKV1ReXq6WlhZlZmZKkgzD0A9/+EOtWLFCGRkZeuaZZ1RaWmrlJYTN0JkCHODw4cPy+/3Mf5sgFotp9+7dhCk4xpAhQ1RXV9e5KADnTu/evWUYhg4dOmR1KZ9Lx8hfbW2totGoampq9JWvfEXBYFCzZ8+Wz+fTXXfdJbfbrXvuuUdXXHGFlixZoq1bt2ratGmqrq7W2rVrFY1G5fF4dPnll8vn82nKlCmqra2Vy+XSX//6Vy1cuLDzzKn3339f9913n6644gq99dZbcrvd+trXvqYf//jHmjRpkpYvX65LLrlEP/vZz2QYhqZOnarLLrtMjz/+uKqqqlhOgS7RmQIcoON+qUR9JzORNDU1qXv37p3vagLJzufzqXfv3qqvr9eQIUOsLiepuVyuzlG/Xr16WV3O59bVYop77rlHFRUVys7O1k9+8hN961vfUn19vY4cOdK5mOK8886Ty+XSqlWrFAgEdNNNN2nRokVKSUnRgAED1NTUpLvuukvDhg1Tbm6uNmzYoN///vfat2+fVq9erZSUFD3zzDO65ppr9KMf/Ug33XSTHn/8cWVmZmrp0qWaN2+efvCDH2ju3Ll6/PHHWU6BLtGZAhyAET/zsHwCTsQSCvMk8hKKkznZYooFCxYoJSVFffv2lc/n00svvaSHHnpI/fv3l9t9/KVrJBKRYRjasGGDfvKTn6hXr17Kzs7WoUOHFAqFdMcdd3R+zEGDBqmhoUHf+ta3dPHFF6u0tFRFRUUKBoPasGGDZs2aJcMwVFtbq+3bt2vOnDnat2+fZs6cqfb2dnXr1q2z3o6PBXQgTAEOwPIJ83C/FJyIMGWeRF2P/lk+vpgiHA5r3bp1ysvLO+E5+fn5WrRokSTp5ZdfVmZmplwul2bPnq0XXnihMxAdPnxYhYWFnR/z6NGjikQievrpp7V582Y99NBDnR9j1qxZ8vl8nX9GU1OTpONdwA0bNsgwDGVnZ5t0FZCICFOAA7AW3Txs8oMTcdaUeZKtM9Xh4yN/I0aM0IwZM3Ts2DHdc889Wr58uerr6zVt2jQdOnRIRUVFevDBBzVmzBjV1dXpvPPO01VXXaXa2lrNnTtXV155pRoaGjo/5muvvaY777xTWVlZ8vv9mjt3rnbt2qX//u//1ty5cyVJffv21d69e/Xyyy+rpKREfr9ft99+u6ZPn676+vrOOuvr6zVw4ECrLhNsiHumAAegM2We6upqTZ061eoyAFMVFBTo73//u9VlOEIyrEfvysyZMzVz5kxJx8f4iouL9fjjj2vgwIGaP3++Fi9erPnz53c+/7HHHtOiRYs0adIkFRcX66qrrtKLL76ojz76SN/4xjd0xx13aMSIEerbt6+qqqqUkpIiSdq9e7dmzZqlbdu2dX6s2bNna9GiRZo3b56OHTum5uZmPfDAA3rttdf06KOP6rrrrtP69euVnZ3N/VI4AZ0pIMnFYjHt3LmTzpRJ6EzBiThryjzDhg1TTU1N0m9P/GSn6tprr+1cTrF8+XJJ0g033HBCp+r++++XdHw5xbXXXquRI0dqxowZeuyxxzqD1Ny5czVp0iRVVFRo0KBBeuqppyRJ8+bN0z//+U8NGzZMq1at0rx58yQdD3gFBQUqKirSTTfdpN/97ncWXA3YmcswjFM9fsoHAdhfbW2tJk+efMKYAs6dQYMGad26dWw1s8C9996rcDise++91+pSHKexsVHjxo3Tvn37rC7FEYYOHarVq1ersLDQ6lIAp+hyHTKdKSDJcb+UeYLBoA4cOKBBgwZZXQpgqn79+unIkSM6duyY1aU4QjKP+gGJhjAFJDnWopuntrZWgwcP7hwnAZzC7XZr6NCh2r17t9WlOEKybvQDEhFhCkhyLJ8wD2dMwclYj26eZN3oByQiwhSQ5BjzMw9hCk5GmDIPY36AfRCmgCTHmJ95ampqOLAXjsVZU+ahMwXYB2EKSGKBQED79u3T0KFDrS7FEehMwcnoTJln8ODBamlpYeEHYAOEKSCJVVVVKT8/n4UIJuGMKTgZZ02Zx+12q6ioiCUUgA0QpoAkxvIJ8xiGoerqasb84Fj5+fmqrq7WZ5xfibOEUT/AHghTQBLjfinzNDc3y+VyKScnx+pSAEt069ZNmZmZHNxrEtajA/ZAmAKSGJv8zNOxfMLl6vKQdCDpsYTCPGz0A+yBMAUkMcb8zMPyCYAlFGZizA+wB8IUkKQMw2DMz0QsnwBYQmGmjjE/7lEDrEWYApLUgQMH5Ha71atXL6tLcQSWTwD/XkKBc69Hjx7KyMjQ3r17rS4FcDTCFJCkuF/KXIz5AYz5mY37pgDrEaaAJMWIn7k6FlAATsYCCnNx3xRgPcIUkKRYPmGeSCSiuro6DRkyxOpSAEsNHjxYTU1NCofDVpfiCKxHB6xHmAKSFGN+5qmvr1efPn3k8/msLgWwVFpamgYOHKja2lqrS3EExvwA6xGmgCTFmJ95GPED/o1RP/Mw5gdYjzAFJKFIJKKamhoVFRVZXYojsHwC+DeWUJinoKBAdXV1CoVCVpcCOBZhCkhCNTU16t+/P2NnJuGMKeDfOGvKPB6PR3l5eYRXwEKEKSAJsXzCXJwxBfwbZ02Zi1E/wFqEKSAJcb+UuRjzA/6NMT9zEaYAa6VaXQCAs6+iokKjR4+2ugzHYAGF9VasWKGHHnpIhmFo/PjxmjNnjtUlORYLKMxVXFys9evXW10G4Fh0poAkxJifOQzDUGFhoQ4cOKAHHnjA6nIcLT8/X4cPH1Zra6sKCwutLsfRHnjgAR0+fFhDhgxRNBq1upykx1lTgLXoTAFJiDE/c7hcLrlcLhmGoSNHjlhdjqONGDFCw4cPV3t7u0aNGmV1OY525MgRGYYhl8ullJQUq8tJeoz5AdaiMwUkmSNHjqi1tVUDBw60uhRHmD59urxerx555BGrS3G8f/3rX3r77betLsPxHnzwQaWnp2vatGlWl+II/fr1UygUUnNzs9WlAI5EmAKSTGVlpYYNGya3m3/eZnjwwQdVUVGhrKwsq0txvJ49e6pXr15Wl+F4mZmZKi8v16OPPmp1KY7gcrlUXFzMqB9gEV5tAUmGET9zeb1eDRkyxOoyAFsZMmQI59yZiFE/wDqEKSDJsHwCAJyFMAVYhzAFJJmKigoVFxdbXQYAwCTFxcWEKcAihCkgyTDmBwDOwnp0wDouwzBO9fgpHwRgL4ZhqFu3bmpoaFB2drbV5QAATHDs2DH17t1bx44dYx09cG64unqAzhSQRBoaGpSVlUWQAgAHycrKUm5urvbs2WN1KYDjEKaAJFFZWalLL71UsVhMq1atsrocAIBJ3njjDUUiEV122WUqLy+3uhzAUQhTQJLIyclRRUWFDh48qG3btlldDgDAJOXl5dq/f78qKyuVk5NjdTmAo3DPFJBE+vfvL4/Ho5qaGg7tBQCHiMViKioq0tGjR3XgwAGrywGSUZf3TKWaWQWAc+uVV17RgAEDCFIA4CBut1tr167lninAAnSmAAAAAKBrbPMDAAAAgLOJMAUAAAAAcSBMAQAAAEAcCFMAAAAAEAfCFAAAAADEgTAFAAAAAHEgTAEAAABAHAhTAAAAABAHwhQAAAAAxIEwBQAAAABxIEwBAAAAQBwIUwAAAAAQB8IUAAAAAMSBMAUAAAAAcSBMAQAAAEAcCFMAAAAAEAfCFAAAAADEgTAFAAAAAHEgTAEAAABAHAhTAAAAABAHwhQAAAAAxIEwBQAAAABxIEwBAAAAQBwIUwAAAAAQB8IUAAAAAMSBMAUAAAAAcSBMAQAAAEAcCFMAAAAAEAfCFAAAAADEgTAFAAAAAHEgTAEAAABAHAhTAAAAABAHwhQAAAAAxIEwBQAAAABxIEwBAAAAQBwIUwAAAAAQB8IUAAAAAMSBMAUAAAAAcSBMAQAAAEAcCFMAAAAAEAfCFAAAAADEgTAFAAAAAHEgTAEAAABAHAhTAAAAABAHwhQAAAAAxIEwBQAAAABxIEwBAAAAQBwIUwAAAAAQB8IUAAAAAMSBMAUAAAAAcSBMAQAAAEAcCFMAAAAAEAfCFAAAAADEgTAFAAAAAHEgTAEAAABAHAhTAAAAABAHwhQAAAAAxIEwBQAAAABxIEwBAAAAQBwIUwAAAAAQB8IUAAAAAMSBMAUAAAAAcSBMAQAAAEAcCFMAAAAAEAfCFAAAAADEgTAFAAAAAHEgTAEAAABAHAhTAAAAABAHwhQAAAAAxIEwBQAAAABxIEwBAAAAQBwIUwAAAAAQB8IUAAAAAMSBMAUAAAAAcSBMAQAAAEAcCFMAAAAAEAfCFAAAAADEgTAFAAAAAHEgTAEAAABAHAhTAAAAABAHwhQAAAAAxIEwBQAAAABxIEwBAAAAQBxSP+NxlylVAAAAAECCoTMFAAAAAHEgTAEAAABAHAhTAAAAABAHwhQAAAAAxIEwBQAAAABxIEwBAAAAQBz+P9kOXeuT95X1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x864 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Stress discretization\n",
    "stress = data[pp.DISCRETIZATION_MATRICES][parameter_keyword][mpsa_class.stress_matrix_key]\n",
    "# Discrete boundary conditions\n",
    "bound_stress = data[pp.DISCRETIZATION_MATRICES][parameter_keyword][mpsa_class.bound_stress_matrix_key]\n",
    "\n",
    "\n",
    "T = stress * u + bound_stress * u_b\n",
    "\n",
    "T2d = np.reshape(T, (g.dim, -1), order='F')\n",
    "u_b2d = np.reshape(u_b, (g.dim, -1), order='F')\n",
    "assert np.allclose(np.abs(u_b2d[bound.is_neu]), np.abs(T2d[bound.is_neu]))\n",
    "\n",
    "T = np.vstack((T2d, np.zeros(g.num_faces)))\n",
    "pp.plot_grid(g, vector_value=T, figsize=(15, 12), alpha=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the traction on face i: T[2*i:2*i+g.dim] is the traction on the face as defined by the normal vectors g.face_normals. This means that for the bottom boundary, the traction T[bot] is the force from to box on the outside (since the normal vectors here are [0,1]), while on the top boundary, the traction T[top] is the force applied to to top faces from the outside (since the normals here point out of the domain)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
